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l. Introduction 
Besa. 

In this monograph, _we givesja mathematical analysis of 
spectral methods for mixed initial-boundary value problems. 
Spectral methods have become increasingly popular in recent 
years, especially since the development of fast transform 
methods 4eee—Gecs—t0)> with applications in numerical weather 
prediction, numerical simulations of turbulent flows, and other 
problems where high accuracy is desired for complicated solutions. | 
We do not discuss the sophisticated applications of spectral 
methods here; a survey of some applications is given in Sec. 15. 
Instead, we concentrate onthe development of a mathematical theory 
that explains why spectral methods work and how well they work. 
Before presenting the theory, we begin by giving some simple 
examples of the kinds of behavior that we wish to explain. 

Spectral methods involve representing the solution to a 
problem as a truncated series of known functions of the inde- 
pendent variables. We shall make this idea precise in Sec. 2, 
but we can illustrate it here by the standard separation of 
variables solution to the mixed initial-boundary value problem 


for the heat equation. 


Example 1.1: Fourier sine series solution of the heat equation. 


Consider the mixed initial-boundary value problem 


2 
watt) = —— (O< x < nm, t > 0) (1. la) 
u(0,t) = u(m,t) = 0 (t > 0) (1. 1b) 
u(x,0) = £(x) (O< x<m). (1.1¢) 


a 


u(x,t) = J a,(t) sinnx, (1.2) 
n=] 
SiAb) ekg A. ARMA Banca (1.3) 
where 
f(x) sin nx dx (tml 2eee¢) (1.4) 


rh 
Ss 
a 
ayn 
ON, 
a 


are the coefficients of the Fourier sine series expansion of 
f(x). Recall that any function in Ly (0,7) has a Fourier sine 
series that converges to it in Lo(0,1); the Fourier sine series 
of any piecewise continuous function f(x) which has bounded 
variation on (0,1) converges to FUE (xt) +£ (x=) ] throughout 


(0,7) (gee Sec. 3). 


A spectral approximation is gotten by simply truncating 


(1.2) to 
N 


Uy (x,t) = J a, (t)sin nx (1.3) 
n=1 


and replacing (1.3) by the evolution equation 


Ga, 2 
ae =e-n an (n=1,...,N) . (1.6) 


with the initial conditions a, (0) i (n=1,.--,N) . 
The spectral approximation (1.5-6) to (1.1) is an ex- 


ceedingly good approximation for any t > 0 as N+ o 


. 


In fact, the error u(x,t) - Uy (% t) goes to zero more rapidly 


2 
“Nt 
than e as N+ for any t > 0. In contrast, a finite 


difference approximation to the heat equation using N grid points 


=2< 


in x but leaving t as a continuous variable (a 'semi- 
discrete’ approximation) leads to errors that decay only 
algebraically with N as N+. [Of course, if we solve 
(1.6) by finite differences in t the error of the spectral 
method would go to zero algebraically with the time step At. 
However, we shall neglect all time differencing errors for now 
and study only the convergence of semi-discrete approximations. 


Time-differencing methods are discussed in Sec. 94] 


Example 1.2: Fourier sine series solution of an inhomogeneous 
heat equation. 

Not all spectral methods work as well as the trivial one | 
just outlined in Example 1.1. Consider for example the solution 


to the problem 


2 
i aa (0< x<m, t> 0) 
ax 


with the same initial and boundary conditions as before. 


The Fourier sine coefficients of the exact solution are now 


e (1.7) 


where e, = 0 if n is even and 2, * 1 if nn is odd. Spectral 


approximations are now given by (1.5) with (1.6) replaced by 


ae li, +~—e (n=1,...,N) ’ 


the solution of which is (1.7) for n#=1,...,N. Now the 
truncation error u(x,t) - Uy (x, t) no longer decays exponentially 
as N+; the error is of order N"> as N*© for fixed 
x, O<x<m, and t > O0O. In other words, the results 

to be anticipated from this spectral method behave asymptotically 
as N*® in the same way as those obtained by a third-order 
finite-difference scheme [in which the error goes to zero like 


Ax? = baaiey Ws For this problem, straightforward solution by 
finite differences may be more efficient and accurate than solution 


by Fourier series. 
The last example may be disturbing but even more serious 


difficulties confront the unwary user of spectral methods, as 
the next example should make amply clear. 

Example 1.3: Fourier sine series solution of the one- 
dimensional wave equation. 

Consider the mixed initial-boundary value problem for the 


one-dimensional wave equation 


n) 


Sulxrt) , Ur exe t (<x <m t 20) — (2,8a) 


u(0,t) = 0 (t > 0) (1.8b) 


u(x,0) = 0 (O< x < 1) (1.8c) 


The exact solution to this well posed problem is u(x,t) = xt. 
This solution can also be found by Fourier sine series expansion 
of u(x,t). To do this, we substitute (1.2) into (1.8) and re- 


expand all terms in sine series. The Fourier expansion of du/dax is 


Qu. J bi (t)sin nx (1. 9) 


ox nad 


where integration by parts gives 


i ie) ae | §u gin nx ax, © - 22 is cos nx dx 
n ¥ ox ™ % 
oo wv 
=~ 2B a(t) f sin mx cos nx dx, (1.10) 
n=l 0 
eo 
oa * 
mtn odd 


+1 
Also the Fourier sine coefficients of x are 2/n(-1)" and 


the Fourier sine coefficients of t are (4t/mn)e,, where e, 0 


if n is even and *,* 1 if n is odd. Equating coefficients 


of sin nx in (1.8a) we obtain 


2 n, 4 - Bit 
<2 (-1)" +9 te, (n a: ree ( ) 


da J nm 


4 
ec. e wa 


m=1 n° -m 
mtn odd 


The Fourier sine coefficients of the exact solution 


u(x,t) = xt are 


a,(t) = - 2 (-1)"t oe ee 


It is easy to verify by direct substitution that these coefficients 


satisfy (1.11) exactly; in particular, the sum in (1.11) converges 


for all t. 


Now suppose we employ a spectral method based on Fourier 
sine series to solve this problem. We seek a solution to (1.8) in 
the form of the truncated sine series (1.4). If the exact co- 
efficients a(t) are used in (1.4) then u(x,t) - Uy (x,t) + 0 
as N-+o ; for each fixed x, 0<x< 1, and t>0O the 
error is of order 1/N as N+ © (see Sec. 3). 

However, it is not reasonable to assume that the expansion 
coefficients a, (t) are known exactly in this case because of 
the complicated couplings between various n in the system 
(1.11, It is more reasonable to determine them by numerical 


sOlution of an approximation to (1.11). Galerkin approximation 


(see Sec. 2) gives the truncated system of equations 


N ‘ 
da 
see t Y a a 2 et + A te, (nai,...,m) (1.12) 
m=l1 n°--m 
m+n odd 


The truncation of the infinite system (1.11) to the finite 
System (1.12) is a standard way to approximate infinite coupled 
systems. Unfortunately, it need not work. In Figs. 1.1-1.2 
we show plots of the approximations Uy (x,t) at t= 5 given 
by (1.4) for N = 50,75. These plots are obtained by numerical 
sOlution of (1.12) with a, (0) = 0; the time steps used in the 


numerical solution of (1.12) are so small that time differencing 


errors are negligible. It is apparent that the approximate solu- 


‘umoys OS[e ST Gai 3B 4x 


=N UOFINTOS Joexe sy_ “eTQTBT[Zeu aie sio11e Buyousiesjtp eweL *(ZT*T) 


JO UOT{eAZSIUT [eX~AewNU Aq pajuye qo sf UOT ANTOS styl “C= 7B OCsN 10; 
(g°T) 03 (a'x)8n uoftjeupxoidde ufzyAeTeD 9y3 jo AOTd y "I'l °3t4 
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Fig. 1.2. A plot of the Galerkin approximation Uy (x,t) to (1.8) 
for N=75 at t=5. This solution is obtained by numerical 
integration of (1.12). Time differencing errors are negligible. 
The exact solution u = xt at t=5 is also shown. 


N=75 a" 


tions with N finite do not converge to the exact solution as N 


increases! The divergence of this spectral method will be ex- 


plained in Sec. 6. 


Not all spectral methods give such poor results. A properly 


ine SY 


formulated and implemented spectral method gives results of 
striking accuracy with efficient use of computer resources. 
4 The choice of an appropriate spectral method is governed by 
two main considerations: 

(i) Accuracy. [In order to be useful a 


Spectral method should be designed to give results 


of greater accuracy than can be obtained by 

more conventional difference methods using similar 
spatial resolution or degrees of freedom. The choice 
‘ of appropriate spectral representation depends on the 


kind of boundary conditions involved in the problem. 


A (ii) Efficiency. In order to be useful the spec- 


q tral method should be as efficient as difference 
methods with comparable numbers of degrees of 
freedom. For similar work, spectral methods 

| should produce more accurate results than 


conventional methods. 


In Sec. 15, we present a catalog of different spectral methods 


Se gies en aR 


and indicate the kinds of problems to which they can be most use- 
fully applied. 


Many examples of efficient and accurate spectral methods will 


be given later. 


scenic smack ne 


2. Spectral Methods 
The problems to be studied here are mixed initial-boundary 


value problems of the form 


Suet) = L(x, tyu(x,t) + f(x,t) (xe D,t20) (2.2) 


0 (x € 9D, t > 0) (2.2) 


B(x)u(x,t) 


g (x) (x e€ D) (2.3) 


u (x, 0) 


where D is a spatial domain with boundary op, L(x,t) isa 

linear (spatial) differential operator and B(x) is a linear 

(time independent) boundary operator. Here we write (2.1-3) 

for a single dependent variable u and a single space coordinate 

x with the understanding that much of the following analy- 

sis generalizes to systems of equations in higher space di~ 

mensions. Also, attention is restricted to problems with 

homogeneous boundary conditions because the solution to any 

problem involving inhomogeneous boundary conditions is the sum of 

an arbitrary function having the imposed boundary values and 

a solution to a problem of the form (2.1-3). The extension to 

nonlinear problems will be indicated at the end of this section. 
Before discussing spectral methods for solution of (2.1-3) let 

us set up the mathematical framework for our later analysis. 

It is assumed that, for each t ,. u(x,t) is an element 

of a Hilbert space 4 with inner product (, ) and norm 

| «||. For each t > 0, the solution’ u(t) belongs to 


the subspace 8B of 4H consisting of all functions ue H 


1 We will often denote u(x,t) by u(t) when discussing u as 
a function of t. 


| 
| 


Rhee” 


satisfying Bu= 0 on 2D. We do not require that u(x,0)=g(x)eB 
but only that u(x,0) eX . The operator I, is typically an 
unbounded differential operator whose domain is dense in, 

but smaller than, pa g For example, if 


L= 0/oax and Me 


Hy 


L5(0,1), the domain of L can be 


chosen as the ‘dense set of all absolutely continuous functions on 
O< x 4 d., 

If the problem (2.1-3) is well posed, the evolution operator 
is a bounded linear operator from ff to B ,  Boundedness implies 
that the domain of the evolution operator can be extended in 


a standard way from the domain of L to the whole space UH 


(Richtmyer & Morton, 1967, p. 34). For notational convenience 
we shall assume henceforth that L is time independent so that 


the evolution operator is exp(Lt). In this case the formal so- 


lution of (2.1-3) is 


u(t) = e&tu(o) +f el (t-8) er cyas (2.4) 
0 


This formal solution is justified under the conditions 
that £(t) , Lf(t) , and L(t) exist and are continuous 
functions of t in the norm ||+*|| for all t20 (see 
Richtmyer & Morton, 1967). 

The seni-discrete approximations to (2.1) to be studied here 


are of the form 


| 
| 
| 


AOE NRT aR ai 


au,, (x, t) 
N mee 
ot N oN 


(x,t) + fy (x,t) (2.5) 


where, for each t , Uy (x,t) belongs to an N~-dimensional sub- 


space 6 N of B , and Ly is a linear operator from H to 6 


of the form 


LL. = P. LP ° (2.6) 


Here Py is a projection operator of Ht onto By and 


fy = Pyf . We shall assume that 8 


N cBy when N<M. 

For definiteness, we shall also assume the initial conditions for 
the approximate equations (2-5) to be uy(0) = Pv (0) where 

u(0) = g(x) is the initial condition (2.3). Specific 

examples of projections Py and the resulting approximations 

Ly will be given below. 

According to this general framework, the formulation of a 
spectral method involves two essential steps: (i) the choice of 
approximation space By and (ii) the choice of the projection 
operator Py - It will turn out that the mathematical analysis 


of the methods also involves two key steps: (i) the analysis of 


how well functions in 4 can be approximated by functions in 


B (see Sec. 3) and, in particular, the estimation of 


N 
{ju - Pu | for arbitrary uel; and (ii) the study of the 
‘stability’ of Ly (see Sec, 4), Finally, there are the 
important practical questions of how to discretize time (see 
Sec. 9) and how to implement spectral methods etficiently (see 
Sec. 10). All these considerations will be reviewed in 


Sec. 15. 


| 


Fa | 


Galerkin approximation 


A Galerkin approximation to (2.1-3) is constructed 
as follows. The approximation Uy is sought in the form of 


the truncated series 
N 
Uy Ore t) = ney in bt) 0, 0) (2.6) 


where the time-independent functions >, are assumed linearly 
independent and One By for all n. Thus, Uy (x,t) necessarily 
satisfies all the boundary conditions. The expansion coefficients 


a(t) are determined by the Galerkin equations 


BE (Opety) | (Oye Uy) + (Oye£) (nel. ee) (2.7) 


or 
da 


N N 
2 (nr Om) ae = J anlo,1bd,) + (o,f) - 
m=1 ea 
These implicit equations for a(t) can be put into the 
standard explicit form (2.4-5) by defining the projection 


operator Py by 


N N 
Pye 2b PrmlOgrt) Oy) (2.8) 


where Pam are the elements of the inverse of the N x N matrix 


whose entries are (died) 


Note that the relation 


N N 


Py = ney ae Pam (me Py) oy, 0) 


holds for all projection operators Py: However, the specific 


projection operator (2.8) is particular to Galerkin approximation. 


The Galerkin equations (2.7) may be characterized as follows 


S. 


. 


al i ta rat os ss 


At each instant t, we assume that the expansion coefficients a_(t) 
n 

in (2.6) are known and seek values for the N independent quantities 

da, /dat (n=1,...,N) that minimize 


WU duy 
(y- baa LU E> = Duy) e 


1 


The resulting equations for da, /at are (2.7). 
Example 2.1: Fourier sine series 


If we choose t = Ly (0,%) and dO) = sin nx , we re-~ 
cover the Galerkin approximations given in Example i.1-2 for the 
heat equation and in Example 1.3 for the wave equation. Every 
function webs (0,1) has a Fourier sine series that converges 


in the L, norm, so that |]u - Pu \]} +0 as N+, 


2 
However, as illustrated by Example 1.3, this does not ensure 


that the Galerkin approximation Uy converges to u aS Nt, 
Example 2.2: Chebyshev series 
LASHYSnev Series 

We choose H to be the space of functions on 
the interval |x| <1 that are square integrable with respect 
to the weight function V/V 1-x2 < If the problem is 


4 +a * £ix,t) (eh < xe 1, b> 0) » (2. 9a): 


t 


u(-1,t) = 0, wlx,0) = g(x) , (2. 9b) 


| which is a sliqht generalization of Example 1.3, it is appro- 
priate to choose the expansion functions for the Galerkin approni- 
mates to be ¢ (x) = 7 (x) - (-1)"T (x). Here T (x) is the 
n n 0 n 
Chebyshev polynomial of degree n definied by T, (cos) = cos no 
2 
when x = cos0; thus, Ty (x) =], T, &) = X, T, (x) # 2x -l, T, (8) 


4x” ~ 3x,... . Observe that OX) satisfies the boundary condition 


6,6) = 0 because Tf) = (-1)" for all n. The properties 


of Chebyshev polynomials are summarized in the Appendix. 


The Galerkin equations (2.7) are obtained explicitly as 
follows. First, the definition of T,, (x) and the substi- 


tution x = cos § imply that 


1 
(TT) = fcosn ® cosm ode= 2°n, 8 
: : 


nm ’ 


; jy 
(£,9) “iy £ (x)g(x)/v 1l-x* dx, 
-1 


Here Cy = 2, Cc, =1 (n> 0) and 6. = 0 if n¢m 1 


ne=m. Therefore, 


T nt+m 
(Onn) =F Say t (1PM ee, 


Next, the Chebyshev polynomials satisfy 


, 
The) ~ T,h- 


ot, ) * <5 


as may be verified by substituting x = cos®. Therefore, 


n(-1)9*4n + 1m n<m, m+n oda 
(¢ 499m) = a (-1)9*1y n>m, m+n odd 


0 n +mM even 


Using these results, (2.7) gives the Galerkin approximation equations 


-15- 


Cea 


da, rae N N 
at et ” me P 
ptn odd 


N 
+ 2(-1)" a 458 % 
(-1) ob Pa, +f, + 2(-1)" £5 (n=l,...,N). 


p odd 


Here f, = (Tf) for n# Q,...,N. 
These Galerkin equations can be simplified by introducing 


N 
the notation a) = - i (-1)™a,e so that (2.6) becomes 
m=1 


N 
uy(xet) = 2 a(t) T(x). (2,10) 
n=0 


Substituting the above expression for Qo» the Galerkin equations 


can be rewritten as 


for a 
n 
da N 
ne 2 ) 4 1 n 
-—> pa. +f + — b(t) (-1)” (n=0,...,N) 2.11) 
at ch p=ntl Pp n c, ¢ oN), ( 
pt+n odd 
N s 2 
J (2 ave o . (2.12) 
n=0 
Here b(t) is a 'boundary' term that ensures maintenance of the 


boundary condition (2.12). Using (2.12) it is easy to show that 


the explicit form of b(t) is 


b(t) = omen 
Kes 


Nv f 2 
[E,crPontagtal| = Sot [as 


n=0 


Tau approximation 


The tau method was invented by Lanczos in 1938 (see Lanczos 
1956). First, the expansion functions o, (n=],2,...) are 
assumed to be elements of a complete set of orthonormal functions. 
The approximate solution Uy (x, t) is assumed to be expanded in 


terms of those functions as in 


N+k 
Y ajlthesi). (2.13) 
n=l] 


Uy (x,t) = 


Here k is the number of independent boundary constraints Bu,=0 


that must be applied. The important difference between (2.13) 


for tau approximation and (2.6) for Galerkin approximation is that 


the expansion functions o, in (2.13) are not required individually 


to satisfy the boundary constraints (2.2). The k boundary 


constraints 


N+k 
ne 4, 56, * 0 (2.14) 


are imposed as part of the conditions determining the expansion co- 


efficients a, of a function in By: Then, the projection operator 


Py is defined by 


where b (m=l1,...,k) are chosen so that the boundary con- 
Mm 


straints are satisfied: BP.u = 0 for alluece}. 

It follows from these definitions that the tau approximation 
to (2.1-2) is given by (2.13) with the k equations (2.14) 
and the N_ equations 


da. 


at = (ob Uy) + (o,f) (n=1,...,N) ° (2.16) 


An equivalent formulation of the tau method is given as 


follows: The equations for the expansion coefficients a, of the 


exact solution u in terms of the complete orthonormal basis ¢, are 


u(x,t) ae & an(tio(x) , 


da, 
= * (¢,/Lu) + (o£) (Hel 2408) (2.17) 


The tau approximation equations for the N+k expansion co- 
efficients of uy in (2.13) are obtained from the first N 
equations (2.17) with u_ replaced by Uy and the k _ boundary 
conditions (2.14). 

The origin of the name 'tau method’ is that the resulting 


approximation Uy is the exact solution to the modified problem 


aty « 
at =L Uy + f + oh Tp ft) onan (Xx) (2.18) 


a i EB ge TT 


which lies in By for all t > 0. For each initial value problem 
and choice of orthonormal basis o, (and associated inner product), 
there is a (normally unique) choice of t-coefficients such that 


Uy € 8, » namely 


TO (Ogos By tt? (p = k+l, k+2, ...) 


The remaining tau coefficients Tyr Toreeer Ty are determined by 
the k boundary constraints 
r) 
to ae 
ot 
and the N dynamical constraints (2.17) for n = 1,...,N. 


Example 2.3: Fourier sine series 


For all of the applications given in Example 2.1, Galerkin 
and tau approximations based on on = = sin nx are identical 
(except for the scaling factor v270 ) since the orthonormal 


expansion functions on satisfy the boundary conditions. 
Example 2.4: Chebyshev series 


If we choose (x) afin T, (x) where c, = 2, 6 #1 


n+l 1m), 0 n 

(n > 0) and apply the tau method to the problem (2.9) the result 
can be recast into the form of equations (2.10-12) with b(t) = 0 
and (2.11) only applied for n = 0,1,...,N-1 instead of 

n = 0,1,...,N. Thus, the tau equations for the one-dimensional 


wave problem (2.9) are (2.10) with 


a N 
=  S wane Pa hy ONES eM oe 
ptn odd 
. n 
= (1) a, (t) = 0 (2.20) 
n=0 


In this problen, & T, (t)= ay - fy while tp (t) = 0 for pri. 


= 
ah lt asl lal aids 


Example 2.5: Laguerre series 


Here we choose ¥ to be the space of functions that are 
square integrable on 0 <x <o@ with respect to the weight 
function e~. We choose the expansion functions to be 
(x) = L(x) where Ly (x) is the (normalized) Laguerre poly- 
nomial of degree an: Lg (x) = |, L, (x) = 1l-x, 

Ly (x) =l-~ 2x + 5 x*, ave 4 


Suppose we wish to solve the problem 


Uy + uy = f(x,t) (0<x<™, t > Q) (2.2la) 


u(0,t) = O, u(x,0) = g(x) (2.21b) 


by seeking an approximate solution of the form 


Uy Ot) = aa a(t) L(x) Cc (2.22) 


To derive the tan equations for a(t) we note that L(x) 
satisfies L, (0) = 1, EL. = 3! = Lae n= 0,l,... and 


n n+l 
- -x > 3 : 
(L-,) 2 Jy L(x) Ly (xe ax Sa . Thus, the tau approximation 


(2.17) is 
da, N 
— © E a, + (Lr £) (n = 0,1,...,N-1l) (2.23) 


p=n+l 


while the boundary condition is 


x 


rm 


Similarly, the Laguerre-tau approximation to the heat 


equation problem 


wy © UY + f (x,t) iG<z2<-s, +t > 6) 
(2.25) 
u(O,t) = 0 u(x,0) = g(x) 
is given by (2.22), (2.24) and 
da. N 
ae (pon-l)a, + (Lyf) (n=0,1,...,Ne1) (2.26) 
p=n+1 


Collocation or pseudospectral approximation 


The projection operator Py for collocation [sometimes 
called the method of selected points (Lanczos 1956) or pseudospectral 
approximation (Orszag 197lc)] is defined as follows. Let 
Xp eXgree ee Xy be N_ points interior to the domain D. These 
points are called the collocation points. Also let o, &) 
(n=1,...,N) be a basis for the approximation space 8 wand 


suppose that det oy, (x) #9. Then for each ue H 


N 
Pu = } an on (®) (2.27) 
n=] 
where the expansion coefficients a, are the solutions of the 


equations 


L ay $(% 4) = u(x 5) ce ee) (2.28) 
n= 


Thus, collocation is characterized by the conditions that 
Pyu(x.) = u(x,) for i= 1,...,N and Piue BN . Notice that 
the results of collocation depend on both the points Xn and 
the functions (x) for n= Ly«ccgN « 


Example 2.6: Fourier sine series 

If we wish to solve the problems formulated in Examples 
1.1-3 by collocation instead of Galerkin or tau methods 
we proceed as follows. We choose the space w= L,(0,7), 
the expansion functions $, (x) = sin nx (n=1,...,N), and the 
collocation points x. = wi/(N+1) (j=1,...,N). The 


J 
collocation equations 


N 
L a, sin ti = (xs) (J=1,...-N) (2.29) 
n= 


have the explicit solution 


N 
2 ' — 
an ~ Nel ja u(x5)sin rs (nels ses oN} — 


This result follows from the relation 


N - 
+. 79M ... tkn _ Nel 
nea sin WoL sin a? 55% 


valid for 0 < j,k < N+l . Thus, 


N 
Pyt = 2 a, sin nx (2.31) 
n=] 


where a is given by (2.30). 


n 
It follows from (2.29-31) that 


N 
PyLPyu = a b, sin nx 
n= 


where b, = - n? an (n=1,...,N) if L= 927ax7 , and 


. wn 
2 m sin yyy 
bsg fae ee (De 
m=1 cos i cos NeL 
mtn odd 


if L= 0/ox. 


Example 2.7: Chebyshev collocation for the wave equation 


Suppose we wish to solve the one-dimensional wave problem 


(2.9) using collocation. Anappropriate basis for the approxi- 


mation space 


F : “2 = tate 
By «is the set of functions o,, (x) Tf) (-1) Ty (x) 


introduced in our discussion of Example 2.2 above. 


(n=1,...,N) 
We choose the collocation points to be the extrema of the 


Chebyshev polynomial T(x) satisfying |x| <1. Since 


Ty (cos @) = cos N6@, these extrema lie at x5 = cos PE for 


j = 0,...,N-l . The point xy = 7 1 is also an extremum of 


Ty («) but it is not included in the set of collocation points 
because the boundary conditions for (2.9) are imposed at x = -] 


so ¢,(-2) =0 forall on. 


As in Example 2.2, the expansion coefficients a, for 
n= 1,...,N may be augmented by defining a=" ‘ (-1)™ an 
b | =) 
: | so that . 
A | PIA N 
Uy(x.t) = 7 ay (tT (x). 
n=0 


It may then be shown that the collocation equations for a(t) 


: that follow from (2.9) are 


N 
at eahad ea, 7c, * e b(t) (-1) (n=0,...,N) (2.32) 
P+tn odd 


N n 
y (-2) a(t) = 0 (2.33) 
n=0 


where f, = (Tf) and & =G =2, col (0<n<N). 


Here b(t) is a ‘boundary’ term that is used to ensure compliance 


with the boundary condition (2.33). It may also be shown that 


1 n 6 N ¢ 
b(t) = - i J. (-1) (nta+é ) = a [sui - } -0)"8] 


The reader should observe the close similarity between the 


Chebyshev Galerkin, tau, and collocation equations for the problem 
(2.9). The only difference between them is the way the boundary 
term b(t) enters. In the Galerkin equations (2.11), b(t) appears 
with the coefficient (-1)"/e.s in the tau equations b(t) enters 
with the coefficient SAN so it appears only in the equation for 
a, as a tau coefficient; with collocation, the coefficient of 

b(t) is (-1)"7e,. This close similarity between the three methods 
for the wave equation can also be seen by observing that when 
f(x,t) is a polynomial of degree N in x, all three approxi- 
mation methods give Nth degree polynomial approximations Uy (x,t) 


that satisfy exactly the initial-boundary value problem 


r) du 
bin + x = f£(x,t) + T(t)Q, (x) (2.34) 


Uy (-1,t) = 0. 


In the tau method, Qhy (x) = Ty (x) 3 in collocation, 


bs =H Er 8 belie + kak 
Qy (x) = Re (x-x5) = 2 i~— T(x) = 2°" (xed) ty (2) 
2 n=0 Sn 


where x, = cos se (4 = 0,...,N-1) are the collocation points; 


finally, the Galerkin equations (2.10) are obtained if 


: (-1)" 
Q(x) = 7 tr (x). 
N net cy n 


BS Os Sr aii idliaes n 
™ EO eee 


For all three methods t(t) is uniquely determined by the 
requirement that Uy (x, t) be a polynomial of degree N in x 


that satisfies the boundary condition Uy (-1,t) = 0.. for allt. 


Example 2.8: Chebyshev spectral methods for the heat equation 


To illustrate further the nature of the differences between 
Galerkin, tau and collocation methods, we apply them to the 


heat equation 
au au, Strep ee Se ee 
ot me . ? 


u(-1,t) = u(l,t) = 0 (t > 0), u(x,0) = g(x) (-1 < x < 1). 
We approximate u(x,t) by 


N 
uy(xrt) = 2 ant) T(x) . 


n=0 


The Galerkin, tau, and collocation equations for a, (t) are all 


of the form 


da N 


tt ce 5; | ae ee 
oo oa p(p°-n")a_+f_ (t)+b, (t)B, +b, (t)B (2.35) 
dt ch p=n+2 pon 1 in ~2 2n 
ptn even 
N N * 
2 & © F ey ae 0, (2.36) 
n=0 n=0 
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where <=. = (Tf). Eqs. (2.36) are just a restatement of 
{) | 
Uy (FL, t) = 0. The terms b, (t) and by (t) in (2,35) are 
boundary terms that ensure compliance with the boundary condi- 
| 
j tions (2.36). The only differences between the three approximation | 
i 
q methods lies in the coefficients B and B,.. 
ln 2n 
In the tau method, 
‘ ei, * Sy N~2 . Bon ‘ SAN i ei 
In the Galerkin method, 
n 
‘<a a Son 
Fig * S ; Bon c ‘ (2.37b) 
n n 
c this result follows using the expansion functions 
Ty (x) n even 
& (x) = P(x} 
4g " Ty) (x) n odd 
( 
that satisfy o, (42) = 0 and augmenting the expansion coef- 
ficients an for n> 2 by <= ) Aon and pes ) Xonel 
\ Finally, with collocation performed at the points x5 = Cos +P 
(j = 1,2,...,N-1) the coefficients B,, and By, in (2.35) 
are given by 
{ 
1 (-1)" ' 
B = = B Se gets = C (2.37c) 
in Sa. 4 2n ¢.. 
It may also be verified that the boundary terms b, (t) 
: and by (t) are of the form 
| 
; 2 N C14 n F 
{ * au a 5.38 
b(t) = cy aM ee a ee ai 4 ato! 1) ty 4 ) 
ax ts okt n=0 a ee 


for i=1,2. Here 


for the collocation method. 


In the previous examples the only difference between Galerkin, 


tau, and collocation approximations is their treatment of the boundary 


terms. ‘However, in more complicated problems, there are significant 


differences between these approximations. The next example illustrates 
the influence of quadratic nonlinearity. 
Example 2.9: Chebyshev approximations to Burgers' equation 


Chebyshev series approximations to the solution u(x,t) to 


Burgers’ equation 


(|x| $1, t>0) 


u(t+l,t) = 
u(x,0) = f (x) 


are obtained by methods very similar to those for linear equations. 


In general, spectral approximations to the nonlinear equation 


(2.40) 


are of the form 


(2.41) 


where Py is a projection operator. The projection operator 


: E : et 
Py can be that for Galerkin, tau, or collocation approximations. 


0 


If we write 


m7 Z 


Uy (x,t) = a(t) PLO), 


then the Galerkin approximation to (2.39) is given by 


we da Pee 2 2 94m 
Cc, <p = -2) , A, ay * * J m(m*-n“)a + b(t) + b_(t) (-1) 
| m| <N m=n+2 
. mtn even 
|p|<n (O<nsN), 
C m+p2n+1 


ntm+p odd 


(2.42a) 


Observe that if (u, Au) = 0 so the system (2.40) has the 
energy integral d(u,u)/dt = 0, then (2.41) has the energy 
integral O(Uyy ¢ Uy) Zt = 0 provided that the projection operator 


Py is self-adjoint. This follows from 


{ (uy PA (Puy) ) = (Paty ACP Uy) ) = 0. 


An example of a self-adjoint projection operator Py is the 


Galerkin operator (2.8). Energy conservation is quaranteed 
only if the inner product used in the definition of the Galerkin 
approximation is the same as that in the energy integral. 


= 29. 


(2.42b) 


where a) = Sim} 2pm} fOr |m| ¢ N. The tau equations are identical 


except that (2.42a) only applies for 0 <n<sNe-2 and by => = 6, 
On the other hand, the collocation equations obtained using the collo- 


Cation points x5 = cos a for j = 1,...,N-l are just (2.42b) and 


m+p>nt+1 M+p>2N-n+1 
ntm+p odd ntmtp odd 


N 
+ vi m(m7-n“) a 
m=n+2 


= = n 
et b, (t) + b_ (t) (-1) (2.43) 


mtn even 
(0SnSN) 


Where €,) = ¢, = 2 and c¢, = 1 for n+ 0,N. Observe the appear- 
ance of the ‘aliasing' term as the second sum on the right side of 


(2.43). Aliasing is discussed in detail by Orszag (197la, 1972). 


Example 2.10: Chebyshev approximations to uy + F(u), = 0 
SURE aR a een , en 
Galerkin and tau approximations to the solution to 


ur + F(u), = 0 (2.44) 


where F(u) is arbitrarily nonlinear, are very unwieldy both 


to write down explicitly and@ to solve on a computer. On the other hand, 
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while the collocation equations may also be hard to write down 


explicitly, they lend themselves to ready solution without their 
eatieik form being known! 
The collocation approximation to (2.44) is obtained as follows. 
We use the relation 
Uy 
(F(uy)), = F* (uy) x (2.45) 
Since du, / ax can be computed explicitly in terms of uy, asa poly- 
nomial in x of degree N-1l, it follows that (F(uy)), can be 
evaluated by this formula at each of the collocation points assuming 
that F'(z) is a known function; thus, the collocation approxima- 
tion to (2.44) is determined. 
There is a slightly different collocation procedure that can also 
be applied to (2.44). It has the operator form 
ee + Py a Py Flug) = 0» (2.46) 
which is usually not the same as the collocation approximation of 
the form (2.41) described above. In this approximation, a hei is 
computed by first using collocation to obtain PF (uy) from Uy amd then 
using the collocation approximation Py/ 9X to 93/3x given in Example 2.7. 
The collocation approximation given by (2.41) ar (2.45) differs from 


(2.46) by the term 
Py Se (I-Py) Flug) 
N 0x N 


which is generally not zero. However, if F'(z) is not known 


accurately then (2.46) may be the only viable method. 


3. Survey of Approximation Theory 


The remarkable convergence properties of spectral methods to 
be discussed later owe to the rapid convergence of expansions of 
smooth functions in series of orthogonal functions. We present 
a summary of the relevant theory here. 

Fourier series 
The complex Fourier series of f(x) defined for 0S x = 2h 


is the periodic function 


pa ikx 
g (x) 7. ae, (3.1) 
where 
an 
aes f t(xye ax . (3.2) 
0 


We shall show below that if f(x) is piecewise continuous and has 


bounded total variation then 


g(x) = $U£(x+)+£(x-)] (3.3) 


for 0<*x <* 2n and g(x) is repeated periodically outside the interval 


0 <x < 2m. In particular, 9(0) = g(2nm) = g(£(O0+)+£(2n-)) . 


The Fourier sine series of a function f(x) defined for 


O<x< nw is the function 


g,(x) = } a, sin kx r (3.4) 


k=] 


where 


= 2 H f(x) sin kx dx. (3.5) 


The Fourier cosine series of a function defined for 0<x<T is 


G(x) = ] ay cos kx , (3.6) 
k=0 
where 
2 1 
a, = To, } f(x) cos kx dx (3.7) 


with cp) = 2, q =1 (k > 0). It follows easily from (3.3) that if 


f(x) is piecewise continuous and has bounded total variation then 


G(x) = £L(x) , (3.8) 
9, (x) = £, (x), (3.9) 


where f(x) = £, (x) = [(f£f(x+)+£(x-)) for O<x<T, 
£ ,(-x) = ~£ (x), £,(-x) = £. (x) for -n <x <0, £ (0) = £,(m) = 0, 
£ (0) = £ (0+), £,(m) = f(n-), and £ (x) and £ (x) are extended 


periodically outside the interval - 1 < x = hs 


Convergence of Fourier series 
To prove (3.3) we define Dy (x) as the partial sum 


9x (x) = } ae a (3.11) 
k=-K 


Using (3.2) and the trigonometric sum formula 


‘ tin sin[ (K+4)s] 
e = 
k=-K sin(¢s) 
we obtain 
1 x sin( (K+3)t] 
Dy (x) - — f (x-t)dt (3.12) 


x-2T  gin[at] 


The kernel sin(Kt#)t/sin4gt of the integral (3.2) is 
plotted for several values of K in Fig. 3.1. This plot 
suggests that when f(x) has bounded total variation the leading 
contribution to the integral as K + © comes from the neighbor- 
hood of t = 0 since the contributions from the rest of the in- 


tegration region should nearly cancel due to the rapid oscillations 
ot the integrand. Thus, 


te sinf (K+e)t] 
£ (x-t)dt (K+) (3.13) 


Dy (x) v SF 


nt sin[st) 


for any fixed e« > 0. Since e may be chosen small we may replace 
sin #t by 4t with a maximum error of ole). Also since f (x-t) 


is piecewise continuous, we may assume that f(x-t) is continuous 


for 05t=c and -e £ t * 0 with at worst a jump discontinuity 


at t = 0. Therefore we may replace f(x-t) by f(x-) for t > 0 


and f(x+) for t < 0 giving 


*0z‘OT 


S =u 203 [3(k4+u)] /(3(4+0)] UTS JO 30Td wv 


\) 
C 
{ 


"t*¢ °*btd 
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ats 


=n edict 


saa eee ND OT 


€ sin(K+34)s 


Gx (x) © [E(xt)+£(x-)] 2 f ds (K+e) | 
0 s 
Since 
€ sin(K+})s (K+dye _. Marke 
2 f —;—_ da = 2 { eons ds % 2 f S425 ds =4 (K+~) 
0 0 Q 
for any fixed ¢ > 0, we obtain 
Gy (x) v lf (x+) +f (x-)] (Keo) + 


proving (3.3). 

In the neighborhood of a point of discontinuity of f(x) 
for x = 0 and x = 20 if £(0+) + £(2m-)] the convergence 
of Gy (x) to its limit (3.3) as K+ is not uniform. To 
investigate the detailed approach of Dy (x) to g(x) near a 
point of discontinuity Xo of f(x), we use the asymptotic 
integral representation (3.13) to obtain 


€ sin[ (K+4)t] 


c 


z 1 
Sy (Xo+ my ~\ 7} 


£ (xg+—— -t)at (K+) 

K+4 
for every fixed z. Since e is assumed small we can approxi- 
mate £ (x 9+s) by £ (xX9+) for O<s<e and by f£(x9-) for 


-e < s < 0. Therefore, for each fixed z and e, 
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We. £ (xp+) z/ (K+) sin (K+4) i £(x9-) € sin (K+i)t 
Ty (Xp oa T ae t c aia sand : 
K : 
z/ (K+) (K+) 
(K+) 
f (x,+) f(x,-) © 
o 2 f sin 5 ds < 2 sin S 4 (K+) 
-e (K+4) ® 
£(x,+) 2 2; £(x,)-) 1 
~ bbe j SB8ax+—O yy Sin as (K+) | 
=—© z 
Since f sin s/s ds = 7, we obtain 
z 1 ; 
Gy (x,t ——) v S[£(x,+)+£(x,-)) + FL£(x,+)-£(x,-)] Si(z) (Kae). 
K‘"0O K+ 0 0 va 0 0 
(3.14a) 
for anv fixed z. Here the sine integral Si(z) is defined 
2 z 
Si(z) = = j sin § as (3.14b) 


A plot of Si(z) is given in Fig. 3.2. 

The result (3.14) shows that if x - Xo = 0(2) as K~+o then 
Sy (x) - Ff£ (xt) +£ (x—-)] = O(1). This shows the nonuniformity of 
convergence of Ty (x) to f(x) in the neighborhood of the discon-~ 
tinuity xX). This nonuniform bemvior of the limit Jx(x)+£(x) as K + * 
is called the Gibbs phenomenon. 

To illustrate the Gibbs phenomenon in an actwal Fourie series, we plot 


in Fig. 3.3 the partial sums of the Fourier sine series expansion 


of the function 


f(x) = x/n (O<x<m) 


The extended function £,(x) is discontinuous at x = 1 leading 


to the Gibbs phenomenon there. 


1.25 


75 


Pe 


«25 


Fig. 3.2. 


A plot of the sine integral Si(z) defined in (3.14b) for 0 < z < 15. 
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As K+, the maximum error of the partial sums of a Fourier 
(complex or sine or cosine) series in the neighborhood of a point 
of discontinuity occurs at the maximum of Si(z) .- 

Since Si'(z) = 0 when 2 = nn for nn ® #1,+2,..., the maximum 
error must occur at one of these points. It is easy to argue that 


the maximum of Si(z) actually occurs at 2#= n where 


> Si(n) # .58949 (3.15) 


Thus the maximum overshoot of the partial sums of the Fourier 


series near a discontinuity occurs near x = XQ + 5 €or & 


K+$ 
large and is of magnitude 
9 (gto) - £ (x9+) v - 08949 [£ (x9+) -£(x9-)] (K+=) (3.16) 
K+ 


where the quantity in square brackets is the jump at Xo: For the 
example plotted in Fig. 3.3 the jump of £, (x) at x = 7 has 
magnitude 2 so the Fourier series gives a local overshoot of 
magnitude 0.179. 

As £+2 0, Si(z) +z 1 so that (3.14) is consistent with 
the convergence of the Fourier series to £ (x9+) just to the right 


of X99 and to £(xX9-) just to the left of x The Gibbs phenomenon 


0° 
Only appears when x + Xo at the rate 1/K as K+», 


Rate of Convergence of Fourier Series 


If f(x) is smooth and periodic, its Fourier series does not 
exhibit the Gibbs phenomenon. The Fourier series of suchan f(x) con- 


verges rapidly and uniformly. Suppose f(x) is periodic and has 


continuous derivatives of order p = 0,l,...,n-l1 and ¢'™) (x) is 
integrable. Applying integration by parts to (3.2), it follows that 
-ikx 


(x)e ax . 


2n 
1 (n) 
a, —_ | f 


2n(ik)™ Jo 


Since £') (x). is integrable, the Riemann-Lebesque lemma 


implies that 


a, << 1/k™ _(k + to). (3.17) 


Note that, because f(x) is periodic, continuity of ¢'P) (x) 
also requires ¢ 'P) (9) = ¢'P) (2m) It follows from (3.17) 
that if f(x) is continuous with £(0) = £(27) and f'(x) is 
integrable then a, << l/k as k+™ ; if, in addition, f£' (x) 
is piecewise continuous and differentiable then a = 0(1/k?) 
as kro. 

Now we can be more precise in our estimates of the error 
Dy (x) = £(x) « If a, goes to zero like 1/k" as k +o 


and no faster, then ¢ (ML) (5) is discontinuous. In this case, 


gy(x) - E(x) = O(4) (K+ @) (3.18) 
K Bhs 
when x is fixed away from a point of discontinuity of ¢ (n-1) 
as K~+o, while 
gy (x) = £(x) = 0(—4y) (K + ©) (3.29) 
K 


when xX - Xp = 0(%) as K+ © where Xo is a point of 
discontinuity of (M21) (9), 

In particular, if f(x) is infinitely differentiable and 
periodic [f(xt+27) = f£(x)] , 9p (x) converges to f(x) more 


rapidly than any finite power of 1/K as K+o@ forall x. 


Fourier sine and cosine series have convergence properties 
very similar to the complex Fourier series just discussed. We 
summarize these properties for Fourier cosine series. If deri- 
vatives of f(x) of order p = 0,1,...,n-l are continuous for 
0 <x < 1 while £'P) (9) = £'P) (¢) = 0 for all odd p<n 
and £ (PD) (x) is integrable, then the Fourier cosine coefficients 


given by (3.7) satisfy 


n 
a, << 1/k Badia I (3.20) 


as may be proven by integration by parts. 


: 
: 
: 

7 


Thus, if f(x) is infinitely differentiable for 0< x<m 
and f!Pt1) (9) = ¢(2Ptl) (7) = 9 for p= 0,1,... then the 
Fourier cosine coefficients ay approach zero more rapidly 
than any power of 1/k as k++ . In other words, if f(x) 
is infinitely differentiable on -~ < x < ~, periodic with period 
2n- (£(x+27) = £(x)) , and even [f£(x) = £(-x)), 
then the remainder after N terms of the Fourier cosine series 
(3.6) goes tO zero more rapidly than any finite power of 1Aj 
as Kw+o , 

To compare the convergence properties of Fourier sine and 
cosine series, we have plotted in Figs. 3.3 and 3.4 some results 
for the Fourier sine and cosine expansions, respectively, of the 
function x/n for O<x<«<T7. As discussed above, the Gibbs 
phenomenon in the sine series expansion is evident at x= 17 (see 
Pig. 3.3). Observe that the error in the N_ term partial sum 
goes to zero like 1/N as N+o when x is fixed 0< x<T7- 
The Gibbs phenomenon near x = 1 slows the convergence of the 
Fourier series for all x. In Fig. 3.4, we plot the error between 
the N term cosine series and x/t . Observe that as N+” the 
error goes to zero like 1/n? for 0< x < nm and like 1/N when 


x = O(1/N) as Ne. 
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Chebyshev polynomial expansions 


The convergence theory of Chebyshev polynomial expansions 


is very similar to that of Fourier cosine series. In fact, if 


g(x) = J a Ty (x) (3.21) 
k=0 


is the Chebyshev series associated with f(x) for -l<x<«l 
then G(@) = g(cos 6) is the Fourier cosine series of i 
F(@) = f(cos 6) for O<6<T. This result follows from 


the definition of TY fx): since T,, (cos 6) = cos n 6, 


eo 
G(@) = g(cos®) = } a,cosné . Thus, 
k=0 


2 i f (cosé) : : ~x?) 73 
= —— cos ké@ dé = os f | £ (x) a, (x) QQ")? ax 
ia F 


(3.22) 


where Cy) =2, CG = 1 (k > QO). 

It follows from: this close relation between Chebyshev 
series and Fourier cosine series that if f(x) is piecewise 
continuous and if f(x) is of bounded total variation for 
~l<x< 1 then g(x) = Sle (xt) +£ (x-) ] for each x (-1 < x < 1) 
and g(l) = f(1-), g(-1) = £(-1+) . also, if £'P) (x) is 
continuous for all |x| < 1 for p = 0,1,...,n-l, and £'") (x) is 


integrable, then 


<< 1/k" (k + ©), (3.23) 
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Since |, (x) | < 1 for |x| <1, it follows that the re- 
mainder after K terms of the Chebyshev series (3.23) is asymptotically 
much smaller than 1/K°") as K+. If f£(x) is in- 
finitely differentiable for |x| <1, the error in the 
Chebyshev series goes to zero more rapidly than any finite 
power of 1/K as K*®@ . 

The most important feature of Chebyshev series is that 
their convergence properties are not affected by the values 
of f(x) or its derivatives at the boundaries x = +1 but 
only by the smoothness of f(x) and its derivatives throughout 
-1< x < 1. In contrast, the Gibbs phenomenon shows that the 
rate of convergence of Fourier series depends on the values of 
f and its derivatives at the boundaries in addition to the 
smoothness of f in the interior of the interval. The 
reason for the absence of a Gibbs phenomenon for the Chebyshev 
series of f(x) and its derivatives at x = +1 is due to the fact 
that F(®) = f(cos 9) satisfies F{2Pt1) (9) = FI4Pt1) (7) 2 
provided only that all derivatives of f(x) of order at most 2ptl 


exist at x= tl. 
An important consequence of the rapid convergence of Chebyshev 
polynomial expansions of smooth functions is that Chebyshev expansions 
may normally be differentiated termwise. Since 


p 
a: ble T(x) = 0 (k2P) (k + ©) 


dxP 


uniformly for |x| < 1 [see Appendix], if a, * 0 faster than any 


finite power of 1/k as k+o then (3.21) implies 
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(as may be proven by an elementary uniform convergence argument). 
While Chebyshev expansions do not exhibit the Gibbs phenomenon 

at the boundaries x = +1, they do exhibit the phenomenon at any 

interior discontinuity of f(x). In Fig. 3.5 we plot the partial 


sums of the Chebyshev expansions of the sign function sgn x: 


T (x) 
n “2n+1 (3.25) 


en wile 


sgn x = 
Near x = 0, a Gibbs phenomenon is observed; for fixed x # 0, 
the error after N terms is of order 1/N. In general, the local 
structure of the partial sums Sy (x) of Chebyshev series near a 


discontinuity of f(x) is, aside from a simple rescaling, given 


| by (3.14): 


z 1 
Gy (X* : )v sl£ (xqt)+£ (XQ-) ] 
[k+3 } i-x? 
+ Af (xpt)-£(xg-)] Siz) (Ke@) 


where |x,| <1 and z is fixed. This equation is derived 
by a simple extension of the argument used to derive (3.14) 


[cf. (3.33) below for the explanation of the origin of the 


scale factor ii=n? Be 
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Rate of convergence of Sturm-Liouville eigenfunction expansions 


Let us consider the expansion of a function f(x) in terms 


of the eigenfunctions 6 of a Sturm-Liouville problem: The 


n 
eigenfunction $,, (x) is a nonzero solution to 


d do, 
Fe P(X get + Cywlx)-a(x)) 6, (x) = 0 (3.26) 


satisfying homogeneous boundary conditions. To be specific in 
our discussion we assume the boundary conditons >, (a) = o, (b) = 0, 
although the analysis applies more generally. We assume that 

< 


p(x) 20, w(x) 2 0, q (x) 20 for a<*x5=b. We will also 


assume that the eigenfunctions are normalized so that they satisfy 


Fwy (8) 8 (8x eo x (3.27) 
and that they form a complete set; the latter property follows if 
An +o as n+ (see Courant & Hilbert, 1953, p- 424). 

The requirement that An +o follows heuristically as follows: 
(3.26) suggests that ¢,(x) has a typical spatial scale of W/V : 
so the requirement that arbitrary f(x) (with arbitrarily small 


spatial scale) be expansible in terms of {>} implies that Ay 


must grow unboundedly with n. 
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We wish to estimate the rate of convergence of the eigen- 
function expansion 
f(x) = } ano, (x) - (3.28) 
n=l 


Using the orthonormality relation (3.27), the L,~- error after 
N terms is 


b N 4s Ce) \ 
f J£(x) - I ao, (x) |? we . a2 ‘ (3.29) 


- n=l 


Thus, the Lj-error may be estimated by calcnlating the rate 


of decrease of an as nwo, 


Orthonormality of {o,} implies that 


b 
a= u £ (x) >, (x) w(x) dx . (3.30) 


Substituting w(x) o) (x) from the Sturm-Liouville equation (3.26) 


gives 
b ad 
aes f (- & Pm BE + a0, ) £O0ax . 
nia 


Integrating twice by parts, we obtain 


b “b 
an = POH) [oy (HD (x)-05 E00) | + Kf Mx o_GwoNax (3.31) 
n 


n 


Nace bait Balen 


= eo df 
nix) = [- Spi $+ qooeoo] /voo. (3.32) 
This integration by parts is justified if f£ is twice differentiable 
and h is square integrable with respect to w. Under these con- 


ditions and recalling that $,, (a) = >, (b) = 0, we obtain 


oh ak ’ ¢ 1 
a,= i [p (a) 9) (a) £ (a) -p(b) 97 (b) £(b) ] + ge 


: . 2 . 2 b 2 
as n+ .®%, since | f how ax | s f h“w dx J ¢° wax = 0(1) as 
a a a. 
ne @ « 


Nonsingular Sturm-Liouville problems 

To proceed further we must distinguish between nonsingu- 
lar and singular Sturm-Liouville problems: a problem is non- 
singular if’ p(x) >0 and w(x) >0O throughout a*x <b. The 
important conclusion from (3.31-32) is that if the Sturm-Liouville 
problem is nonsingular and if f(a) or £(b) is nonzero then 


7 oh = [pla)g,(a)f(a)-plb)gn(b)£(b)] = (n + @) 


' 
(Notice that if $,, (a) = 0, then $,, (x) = 0 since (3.26) is 


second-order differential equation and p(x) # 0). 


It is well known [Courant & Hilbert 1953] that the asymptotic 
behavior of the eigenvalues and eigenfunctions of a nonsingular 


Sturm-Liouville problem are given by 


aT 


, [Ps vE ax]? (n + ©) (3.34) 


x 
6, (x) vA, sin(yh,, J VE ax) (n + @) (3.35) 


Using these relations in (3.33), we find that a, behaves like 


1 
_ > ©, 
n as n 


This behavior of an leads to the Gibbs phenomenon in the 
expansion (3.28) near those boundary points at which f(a) or 
4 f(b) # 0. The asymptotic behavior (3.34-35) implies that this 
q Gibbs phenomenon is asymptotically identical to that exhibited by 


Fourier sine series provided we use the stretched independent variable 


b 
X = n(x-a)¥wla)/pla) /{f /wls)7pl(s) ds (3.36) 
a 


1} near x = a anda similarly stretched coordinate near x = b. 
If f(a) = f(b) = 0,then a, << l/n as n+, However, a 
further integration by parts in (3.31) shows that if the Sturm- 


Liouville problem is nonsingular and if h(a) or h(b) #0, 


then a, behaves like i, as n-+o, In general, unless. f(x) 


ee 


satisfies an infinite number of very special conditions at x =a 
and x = b, then ay, decays alqebraically as n>. 
These results on algebraic decay of errors in expansions 
based on nonsingular second-order eigenvalue problems generalize 
to higher-order eigenvalue problems. For example, as n+~, the expansion 
o 
coefficients in a, in f(x) = abo ao, (x), where {$ (x)} are 


the normalized 'beam' functions 


OAT = Oy On(t1) = ON (tl) =O, 
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behave like + if £(#1) # 0 (implying a Gibbs phenomenon at the 


boundaries x = al); Like 4, if £(t1) = 0 but £'(#1) #0, 
n 
like + if £(t1) = £'(41) = 0 but £''''(+1) #0, and so on. 
n 


Singular Sturm-Liouville problems 
If p(a) = 0 in (3.33) then it is not necessary to require that 
' 
f(a) = 0 to achieve Sh ge (ERR For this reason, 


n 
expansions based on eigenfunctions of a Sturm-Liouville problem that 


is singular at x =a do not normally exhibit the Gibbs phenomenon 
at x =a. Furthermore, if the argument that led to (3.33) can be 
repeated on h(x) given by (3.32) [this is possible if p/w, p'/w, 
and g/w are bounded and all derivatives of f are square integrable 
with respect to w] then the boundary contribution to an from 
x = a is smaller than ri as n+o, If there are also no 
boundary contributions ss from x = b_ when the operations leading 
to (3.33) are repeated indefinitely [which is true if p(b) = 0], 
then an decreases more rapidly than any power of 7 as n+o, 
The important conclusion is that the rate of sccieyenss of 


eigenfunction expansions based on Sturm-Liouville problems that are 
Singular at x = a and at x = b_ converge at a rate governed by 
the smoothness of the function being expanded not by any special 
boundary conditions satisfied by the function. 


Fourier-Bessel series 
A Fourier-Bessel series of order 0 is obtained by choosing 
the expansion functions to be the eigenfunctions of the singular 


Sturm-Liouville problem 


ad, Mn 
ax* ax + AnXon = 0 (0 < x < 1) (3.37) 


,(1) = 0, 6. (0) finite. 
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Therefore, p(x) = w(x) =x in (3.26) so the problem is singular 


at x= 0, but nonsingular at x = 1. The eigenfunctions are 


o,(x) = Ig (Ign%) 


where Jo is the Bessel function of order 0 and Jon is its 
2 
nth zero, Jo (Igy) = 0. The eigenvalues An = Jon satisfy 


Sai ~ (n- 4) (n+0) . 


The Fourier-Bessel expansion of a function f(x) is given by 


x). (3. 38a) 


ao 
£ (x) “2 Qn JQ (Jon 


The expansion coefficients a, are given by (3.30): 


1 
2 
a, = ——+ | tt(t)og(igtiat. (3. 38b) 
59 Gon) 


because 


fea (j.t)2at = 4a%(5_)? 
0 0 ‘“on 0 *-on : 
For example, the Fourier-Bessel expansion of f(x) = 1 is 


1=- 


J,(3_ x) (3.39) 
’ 0 ‘“on 
n= Jon70 Son) 
In Fig. 3.6 we plot the 10, 20, and 40 term partial sums of the 
series (3.39). There are three noteworthy features of this 
plot: 
(i) At x = 1 there is apparently a Gibbs phenomenon. In 


fact, it is easy to show that this Gibbs phenomenon has the same 


structure as that for Fourier sine series: 
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Fig. 3.6 A plot of the partial sums of the 
Fourter-Bessel series expansion (3.36) of the 
function 1 truncated after N=10,20,40 terms. 


Observe the Gibbs phenomenon near x = 1. Also 


observe that the series converges like 1/N for 
fixed x satisfying 0<x< 7 and like 1//N 


for x near 0. 


| + 
\This behavior is not too surprising because JQ (z)v(2/n2) cos (z~41) 


as z++~, so’the large n_ behavior of (3.39) can be asymptotically 


approximated by that of Fourier series. 


(ii) For fixed x satisfying 0< x <1, the error after 
N+l terms of (3.39) is 
N 2 


1 
1+ Ty (Gon®) = 0%) 


n=O J6n99 on) 

In fact, the nth term of (3.39) has magnitude of order l/n 

and oscillates in sign roughly every min (, rc) terms. The 

error in such an oscillating series is ‘of order 1/N after N terms. 
(iii) At x = 0, the series converges (so there is no Gibbs 

phenomenon there) but the convergence is very slow and oscillatory. 

In fact, the error after N terms is of order (21) St1 vq ° 


This follows because 


vv ; (-1)" rn (-2) N*2 
n=N+l “nn ¥ 2N 


(3.40) 


This slow rate of convergence near x = 0 holds even though the 
eigenvalue problem is singular at x= 0. There are two reasons 

for the slow convergence of Fourier-Bessel series near x=0. First, 
the Gibbs phenomenon at x = 1 affects the rate of convergence 
throughout 0 *x £1. In fact, this is the sole source of the 


‘behavior (3.40). However, when f£'(x) # 0, slow convergence near 


x = 0 can also originate because p(x) = w(x) = x gives 
p'/w = 1/x which is singular at x= 0 so h(x) given by (3.32) 
is singular at x=0 if £'(0) +0. 


Chebyshev series revisited 


Chebyshev polynomials are the eigenfunctions of the singu- 


lar Sturm-Liouville problem (3.26) with p(x) = fine’, 
=x? 


w(x) = 1/vl-x" , q(x) = 0, - 1 S$ x S$ 1, and the boundary conditions 


that (+2) be finite. The eigenvalue corresponding to TW *) is 


An = n*. Since p/w = inx* and p'/w = -x are both finite for 
|x|si, it follows that the argument leading from (3.30) to (3.33) 
can be repeated on h(x) given by (3.32) so long as f(x) is 
sufficiently differentiable. Therefore, the Chebyshev series 
expansion of an infinitely differentiable function converges 
faster than any power of 1/n as n+, as shown fdllowtng (3.23) 
by a different method. 


To illustrate the convergence properties of Chebyshev series 


expansions, we study the rate of convergence of the series 


sinM-(xta) = 2) 2, (Mn) sin (Mratann)?, (x) ix] £2 


n=0 “n 

(3.41) 
Since J,, (MT) + 0 exponentially fast as n increases beyond Mn, it 
follows that (3.41) starts converging very rapidly when more 
than M™ terms are included (see Fig. 3.7). This result leads to 
an heuristic rule for the resolution requirements of Chebyshev 
expansions. Since sin Mm(x+ta) has M complete wavelengths lying 
within the interval |x| < 1, we argue that Chebyshev expansions 
converge rapidly when at least m polynomials are retained per 
wavelength. In general, we expect that the Chebyshev expansion 
of a function that oscillates over a distance \ converges rapidly 


if 27/X polynomials are retained. Fewer polynomials are required 


only (see below). 
if the region of rapid change of the function occurs,at the boundary, 
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Fig. 3.7. A plot of the Lo-error in the Chebyshev series expansion (3.38) of 


sin(Mmx) truncated after Ty @) versus N/M. The various symbols represent: 
OM= 10; <M= 20; AM= 30; OM = 40. Observe that the L,-error approaches 


zero rapidly when N/M > 7. 
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The Chebyshev polynomial expansion of a function f(z) 


that is analytic in a region of the complex-z plane that includes 


the interval -1 < z S$ 1 converges at least exponentially fast 


as nro, If f(z) has singularities in the finite-z plane 
then 


F 1/k 1 
lim sup la = = 3.42 


where R is the sum of the semi-major and semi-minor axes of 


the largest ellipse with foci at z= +1 within which f(z) 


has no singularities. Thus, the Lj-error (3.29) after N terms of 


the Chebyshev expansions decays to 0 roughly like ih as N7™, 


To prove (3.42), we note that 


2 ae f(z) T, (z) 


a = dz 
n Te, ay fy 22 
-1/2 -n 
» —<- [ e£(2) (1-24) (2+ ¢z2-1) az (3.43) 
te, Je 


where C is any contour that encircles the interval (-1,1) 

just once and does not enclose singularities of f(z). 

Eq. (3.43) follows because 2T, (2) = ({z + fz =i)" 4 (z - vz ee } 
where we choose the branch of Jo2-1 satisfying for. Vv Z 

as z2+@. Since (2 + /z*-1)"+ 0 as z+ with this 
choice of branch cut, we can expand the contour C to infinity 

by Cauchy's theorem and pick up the contributions from the 
singularities of f(z). If the 'nearest' singularity is a pole 
at z= Zo with residue r (other singularities may be treated 
similarly), then 


. = r (29 + fy? = 1)7* (n+), 


To complete the justification of (3.42) we need only show 


that [Zo + v22-1| = R. Recall that an ellipse with foci 


at +1 satisfies x7 /a* + y7/B" = 1 with a2 - B? =l. 


If 2z lies on this ellipse, then setting 2, = A cos 6 + iB sin 6 


0 0 
it follows that zo + z= = (A+B) e*® = Re*? r 

Let us give an example of the behavior (3.42). The function 
f(z) = tanh (10 z) has poles at z= +in/20 . Thus, 
R = 7/20 + vl + (7/20) 2 = 1.16934. |The Chebyshev expansion 
cOefficients of f(z) satisfy aon = 0 (because f(z) is an 


odd function), while a, # 1.2679, a, + - 0.4089, a, = 0.2300, 
and so on. The rms (L,) error ey [see (3.29) ] 

obtained by truncating the series for f(z) after Ty (2) 
satisfies (e,/e,,) = (1.175) 2 > C47/eg9 + (1.16935) 2 ’ 
demonstrating (3.42) for this case. The error en is smaller 
than 0.01 for N > 25, which again illustrates the result that 
roughly wm polynomials per 'wavelength' are required to resolve 
a function; the function f(z) has a region of rapid change 
near x = 0 of width roughly 0.1. 

If f(z) is entire, R= in (3.42) so its Chebyshev 
expansion coefficients decay faster than exponentially. More 
precisely, the method of steepest descents applied to (3.43) gives 
the following result: if f(z) is entire and 


£(z) = 0(|z|8 exp \z|*) as 2+, then 
lim sup(tnja,|)/(nénn) = - L (3.44) 
no 


For example, sin Mn(zta) is entire with a= 1 while its 
Chebyshev coefficients in (3.41) satisfy an = o( (Mn) "/nt) 
as n+, in agreement with (3.44). Also, a polynonial has 


Chebyshev coefficients that satisfy (3.44) with a=0. 
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Finally, we remark that the Chebyshev series expansion 
(3.21-22) of an arbitrary function g(x) has a maximum 
pointwise error that does not differ drastically from the 
smallest possible maximum pointwise error of any Nth degree 
polynomial, the so-called minimax error. In fact, the maximum 
pointwise error of the Chebyshev series (3.21) truncated after 


2 gnN) times larger than the minimax 


TY (x) is at most 4(1 +7 
2 


error (Rivlin 199). Since 4(1 +7 “ &nN) < 10 for 
N < 2,688,000, the Chebyshev series is within a decimal place 
of the minimax approximation for all such polynomial approxima- 


tions. 


Legendre series 


Legendre polynomials are the eigenfunctions of the singular 


Sturm-Liouville problem (3.26) with p(x) = jan", q(x) = 0, 


w(x) = 1 for -1 Sx <1 and the boundary conditions are 


An = n(n+l) and its eigenfunction is $,, (x) = P(X) the 

: 2 
Legendre polynomial of degree n. Since p/w=1-x and 
p'/w = -2x are both finite for |x| <1, it follows that the 


Legendre series expansion of infinitely differentiable functions 


converges faster than algebraically. 


To illustrate the convergence properties of Legendre series, 


we study the convergence of the series 


sinMn(x+a) = — wag 2h neg Hm) sin (Nra+ant)P. (x) 
(3.45) 


Since the expansion coefficients in (3.45) vanish rapidly as n 


increases beyond Mt, we conclude that Legendre polynomial expansions 


r- 


o£ smooth functions converge rapidly provided that at least 


t polynomials are retained per wavelength, (see Fig. 3.8). 


When a discontinuous function is expanded in Legendre series, 
the rate of convergence is no longer faster than algebraic. In 
the neighborhood of a discontinuity, a Gibbs phenomenon occurs 
whose local structure is the same as that for Fourier series 
with a suitable stretching of the coordinate. For example, the 


Legendre series expansion of the sign function sgnx is 


sgnx = ) (-1)" (4n+3) (2n)! 


P (x) (3.46) 
n=0 22721 (nga) in! 2n+1 


The partial sums of this series are plotted in Fig. 3.9. Three 
features are noteworthy: 

(i) The Gibbs phenomenon near xX = 0 has the same structure 
as that for Fourier series. 

(ii) The error after N terms behaves like 1/N for |x| <1, 
x #0. This follows from the fact that the (2n+1l)st Legendre 


coefficient in (3.46) satisfies 


(4n+3) (2n)! 1 
a = (-2)7 —eneeiisn) = 0(—) (n+0) (3.47) 
n 22Mt1 (n41) int “nr 


and the estimate 


P(x) = 0(—) (n+ ©) 
vn 


for |x| <1; the series (3.46) is analtemating series if x 


is fixed away from zero so the error after N terms is at most 


of order a_P =o (-2) * é 
nn n 
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Fig. 3.8. A plot of the Lj-error in the Legendre series expansion (3.39) of 
sin(Mmx) truncated after Py) versus N/M. The various symbols represent: 
O M= 10; <*M= 20; AM= 30; OM = 40. Observe that the L-error approaches 


zero rapidly when N/M > 7. 
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(iii) The series converges only like 1/YN at x= #1. This 


follows from (3.47) because P(t) = (+1)” for all n. Thus, 


an interior Gibbs phenomenon in a Legendre series expansion has 


a ‘long-range' effect in the sense that it seriously affects the 


rate of convergence at the endpoints x = +l of the interval. 


In contrast, the error of the Chebyshev expansion of sgn x 


plotted in Fig. 3.5 decay like 1/N at x= #1. This behavior 


is quite general; the boundary errors of Legendre polynomial 


expansions decay to zero roughly a factor YN slower than the 


boundary errors of Chebyshev expansions. 

The rate of convergence of Legendre expansions of a general 
function f(x) may be estimated as for Chebyshev expansions. 
In particular, the results (3.42) and (3.44) hold provided that 


f(x) satisfies the stated conditions and (3.23) holds with 


only minor modifications. 


Resolution of thin boundary layers 


Legendre and Chebyshev polynomial expansions give an 


exceedingly good representation of functions that undergo rapid 
changes in narrow boundary layers. Consider the sequence of 
functions Gs (x) = f(x) exp[(x-1)/5] as §6+0 with Reé>0 for 
a fixed smooth function f(x). As 6>0, gs (x) develops a 
boundary layer of width § near x=l. It may easily be shown 


that the Chebyshev expansion coefficients of Tg (x) satisfy 
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12 
a, ~ (26/n)1/? £(1)e°2 ™ § (noe 5 bn*=0(1)) (3.48) 


provided that Re6é>Q, Thus, if N_ polynomials are retained, 


the rms error e¢ in the Chebyshev expansion of Gs (x) satisfies 
1 2 
mMe~- x (Reg) N (N+), (3.49) 


The result (3.49) implies that as 6+0, the number of 
polynomials required to reach a specified error bound increases 
only as 1//§, in contrast to a uniform grid representation of 

5 (x) that would require order 1/§ grid points in the interval 
|x|<1. In fact, to achieve 1% maximum pointwise error in boundary 
layers of thickness § at the ends of the interval -I<x<l, it 


is necessary to retain only 


N ~ 3//YRe 6 (3.50) 


polynomials as §>0. 


Heuristically, the reason that Chebyshev expansions represent 


boundary layers so well is that the extrema of T, fx) occur 
2 2 
= cosnj/n for j70,1,.../n. Since xX g-x,~ 7 /2n and 


. n+ it follows that these polynomials can 


at *5 
-X_ ~ n?/2n as ’ 
4 -2 


x 
n=1 
resolve changes over distances of order n 


The convergence properties of Legendre polynomial expansions 
of boundary-layer functions are similar to those of Chebyshev 
expansions. In particular, (3.49) and (3.50) are both still 
valid. In Fig. 3.10 we compare the spatial distribution of the 
errors in Chebyshev and Legendre polynomial expansions of the 


function g(x) = el 0 (x-1) | 


which has a narrow boundary layer 
of width 1/100 near x=l. Apparently for x away from the 
boundaries x=+l, the Legendre expansion has somewhat smaller 
errors, while near x=+l the Chebyshev expansion has smaller 
errors. 


The Legendre expansion gives the polynomial Qyy (*) of 


degree N_ that minimizes 
1 
2 
Sf tatxy - Qsy (0) | dx 
1 


while the Chebyshev expansion gives that Qy that minimizes 


1 
f [old = Q. 00 |? (rex?) 727? ax, 


The Chebyshev expansion also gives a smaller maximum error 


max |g(x) - Qyy (0) | 
|x| <1 
than the Legendre expansion by roughly a factor 2/vN ; 
as remarked above, the Chebyshev Qhy (2) is usually remarkably 


close to the minimax polynomial that minimizes the maximum 


error. 
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Laguerre polynomials are the eigenfunctions of (3.26) with 
p(x) = xe™*, q(x) = 0, w(x) =e* for 0< x < @ with 


1 
— th 


e $,, (x) bounded at x = 0 and », The n~ eigenvalue is | 
a. =n and the associated eigenfunction is $9, (x) = L(x), 
the Laguerre polynomial of degree n. If f(x) and all its 


derivatives are smooth and satisfy 


f(x) = O(e°*) (x + ©) 


for some a < 4, it is easy to show by retracing the derivation 


of (3.33) from (3.30) that the Legendre expansion 


f(x) = } a, L(x) 
n=0 


converges faster than algebraically as the number of terms N > ©. 


To illustrate the rate of convergence of Laguerre series, 


we consider the expansion of sinx: 


sinx = - sw72 coal Fins jf, 0 (3.51) 


n 


which converges for all x, 0 x < ©, Since 


L(x) % sah ef X74, -t cos(2Ynx -in] , 
VT 


{see Erdelyi et al 1953, Vol. II, pg. 200) it follows that if 


N > > x, then the error after N terms at x is roughly 


etx 


2N72 (gx) 2 


This error is small only if N&n2 >x or N ~*~ 1.44x. Since 
the wavelength of sinx is 21 , Laguerre expansions require 


| approximately 9.06 polynomials per wavelength to achieve high 
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accuracy. (This figure may be reduced to about 6.53 polynomials 
per wavelength by using the modified Laguerre expansion 
) aL, (x)eo* and optimizing the choice of a.) Thus, Laguerre 
expansions require many more terms to resolve a function of given 
complexity than do either Chebyshev or Legendre expansions. The 
reason is that significant weight is given to x ++o© in the 
Laguerre series where sinx has an essential singularity. 

In Figs. 3.11-13, we plot the partial sums of (3.51) with 
N = 10, 20, and 40 terms. Observe that the number of wavelengths 


of sinx represented accurately by (3.51) is roughly N/9. 


Hermite expansions ; 
Hermite polynomials satisfy (3.26) with p= 63% q(x) = 0, 
2 


] 
| 


2 
w(x) =e * for -°2<x<», $, (x)e #x” pounded as |x| + ©. 


The Hermite polynomial H,, (x) of degree n is associated with 


dh. = 2n. If £(x) and all its derivatives satisfy 


the eigenvalue = 


2 
f(x) = 0(e%™ ) (|x| + ~) 


for some a < 4, then the Hermite expansion 


f(x) = ) an H, (x) 
n=0 


converges faster than algebraically as the number of terms N+ ~. 


This is proved by retracing the steps leading from (3.30) to (3.33). 


To study the rate of convergence of Hermite series, we consider 


a ee ee ET TIT aT EE 


the expansion of sinx: 


° 1 
sin x = — H (x) (3.52) 
} 22ntt onery, = 2ntt 


n=0 
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Fig. 3.11. A plot of sin x and the partial sum of the 
sin x in (3.5) truncated 


Laguerre series expansion of 
after the Laguerre polynomial Ly ®) with N=10. 

Observe that the accuracy of the approximation decreases 
as x increases and that roughly one wavelength of sin x 
is approximated accurately by this truncation of (3.51). 


*x ufZs jo 


syyZueTeaRPa Z ATYBSno1 10zZ 93eINDOe St uoTIeW;xXOAdde ey3 3eY3 
2A198q0 °OZ=N UITA (x) Sq 103je pezeouni, St ({>¢°€) seties 
aliensey] sy2 Wey? Ydeoxe T[°¢ “3TqZ se owes “ZTE °3TA 


N=40 


“a, 
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Fig. 3.13. Same as Fig. 3,12 except that the 

Laguerre series (3.51) is truncated after Ly) 
with N=40. Cbserve that the approximation is 
accurate for roughly 44 wavelengths of sin x. 
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Since the asymptotic behavior of H,, (x) is given by [Erdelyi, 


et. al 1953, vol. II, pg. 201) 


2 
H,, (x) \ e?* aT cos (/Y2ntl x - énmT) 


as n-+o for x fixed, it follows that the error after : 
terms of (3.52) goes to zero rapidly at x only if N 2 wx . 
This result is very bad; to resolve M wavelengths of sinx 
requires nearly mu? Hermite polynomials! [By expanding in the 
series /J a, H,(x)e °* -and optimizing the choice of a, it is 
possible to reduce the number of required Hermite polynomials to 
about 31 = 7.85 per wavelength, but this is still quite poor.] 
Because of the poor resolution properties of Laguerre and 
Hermite polynomials the authors doubt they will be of much prac- 


tical value in applications of spectral methods. 


4. Review of Convergence Theory 


The fundamental problem of the numerical analysis of 
initial value problems is to find conditions under which 
Uy (x,t) converges to u(x,t) as N+ © for some time in- 
terval 0 < t<T and to estimate the error |lu - Uy || ; 
The principal result is the Lax-Richtmyer equivalence theorem 
which states that stability is equivalent to convergence for 
consistent approximations to well-posed linear problems. The 
terms stable, convergent, and consistent relate to technical 
properties of the approximation scheme which are defined below. 


An approximation scheme (2.5-6) is stable if 


Lyt 
lle if <° Ret) (4.1) 


for all N where K(t) is a finite function of t. Here 


the operator norm is defined by 


t N 
We | = max tle "ull 
uew lull 


An approximation scheme is convergent if 
jJu(t) - uy, (e) II + 0 as N+ © 


for all t inthe interval 0<t<T and all u(0)e} and 


f(t) eX. Finally, an approximation scheme is consistent if 


[uu - yall + 0 
(4.2) 


_ } > 
Ju - Pit * 0 


as N+o for all u ina dense subspace of Af . 


The classical Lax-Richtmyer equivalence theorem relating 
the above definition states that "a consistent approximation to 
a well-posed linear problem is stable if and only if it is 
convergent." In this monograph we are confronted with some 
subtleties regarding the notions of stability and convergence. 
Because a precise understanding of the ideas of stability and 
convergence is important to the theory of algebraic stability 
given in Sec. 5, we outline here the proof of the equivalence 
theorem. 

Proof of the Equivalence Theorem 
To show that stability implies convergence we use (2.1) and 


(2.5) to obtain 


Au-u,, 
oe ol Ly (u-uy) + Lu - LY oe fy : 
Thus, 


Lyt 
u(t) - uy (t) =e [u(0)-u, (0) ] 


t L,(t-s) 
ait Pik [Lu(s)-Lyu(s)+£(s)-£,(s)] ds, (4.3) 
0 


Using (4.1) and (4.3) and the triangle inequality we obtain the 
estimate 
Ju(t)-uy(t) || < K(t) Ju (0) =u, (0) || 


t 
+ f K(t-s) |||Lu(s)-L,u(s) || + Il £(s)-£, (s) |I ds 


(4.4) 


Thus, if u(t) belongs to the dense subspace of }€ satisfying 
(4.2) and if f(t) belongs to the dense subspace of )€ satisfy- 


ing ||f - Pyfil +0 as N+, then Jat) - uy(t) || + 0 
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as Ne. Since all solutions u(t) of (2.1) can be 
approximated arbitrarily well by functions satisfying (4.2), 
the proof that stability implies convergence is completed. 
Conversely, to show that convergence implies stability, 
we first observe that, for any ue, lle Na is bounded 
for all N and each fixed t. In fact, convergence implies 


Lt Lit 
o < |lle % uff - ffleb ull] < lleXu-eMtul] + 0, (N>@) 


while well-posedness requires that ertul is finite. How- 
ever, waxite ™ ui may depend on u andon t,_ so stability 
is not wad proved. To complete the proof we use the fact that 
a is a Hilbert space. The principle of uniform boundedness 
(Richtmyer & Morton 1967) implies that if le § ull is bounded 


Lt 
as N+ for each t and ue then |le % || is bounded as 


N+ oe for each t. This proves stability and completes the 
proof of the equivalence theorem. 

Using the equivalence theorem, the study of the convergence 
of discrete approximations to the solutions of initial-value problems 
is reduced to the study of the stability of the discrete approxima- 
tions, assuming the approximations are consistent. Thus, the de- 
velopment of conditions for the stability of families of finite- 
dimensional operators Ly is of primary interest in numerical 
analysis. 
Von Neumann Stability Condition 

The simplest condition for stability is due to von Neumann. 
Let us suppose that the Hilbert space +’ possesses the inner product 
ie Using the inner product, we define (neglecting the compli-~ 


* 
cations due to boundary conditions) the adjoint L of an operator 


L as that linear operator that satisfies 
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E | sional approximation Lye the matrix representation of L is 
the adjoint of the matrix representation of Ly (see Sec. 2). 


(u,Lv) = (L u,v) for all u,v in }¢ . For the finite dimen- | 
The operator Ly is said to be a normal operator if Ly commutes 
i 
| 


bd] bi * 
with Ty 80 Uyly = Ely - 
The von Neumann stability condition is that stability of 


normal operators Ly is equivalent to the condition 
Re Ay << CC (4.7) 


where AN is any of the eigenvalues of any of the operators 

Ly and C is a finite constant independent of N. To prove 
* 

this, we note that if L is normal, then Ly and Ly as 


N 
well as exp (Lt) and exp(Lyt) are simultaneously diag- 
nolizable. Therefore, 


2(ReA,)t 


= max e 


tt 
Nv), 2 
N 


= max 
ueH 
where AN are the eigenvalues of Ly ‘ Thus, the von Neumann 
condition (4.7) is equivalent to the stability definition (4.1) 
with K(t) = exp(2Ct) . 
The von Neumann condition gives an operational technique 
for checking stability of normal approximations: compute the 


eigenvalues of Ly and check that the real parts of the eigen- 


values are bounded from above. 


Example 4.1: Symmetric hyperbolic system with periodic 
boundary conditions 


Let us apply the theory just discussed to the stability 
of difference approximations to the m-component symmetric 


hyperbolic system 


du(x,t) . , 2ulx,t) 
— A = (4.8) 
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ey 


SaaS el 


: . > 
with periodic boundary conditions a(0,t) = u(l,t) - 
Here u is an m-component vector and A is a symmetric 
mx m matrix. 


If we discretize in space using second-order centered 


differences, we obtain 


du. u - u. 
eee: j+1 j-1 ae 
at A Thx Gj Weg2 7s cceteg iu) (4.9) 


Uy (t) = wy fe) ’ u, (t) = Uy 4, ft) 


where u, (t) = u(k/N,t) and Ax =1/N. The system (4.9) 
is equivalent to the system of mN equations 


yg Cn weeny 
+: Bou (4.10a) 


where a is the column vector whose transpose is 


“a 


uP = (8) ,0,,..-i)). Here B is the mN x mN matrix given 


by the Kronecker product 
Be A@obdD =, (4.10b) 


where A is the mxm matrix in (4.8) and D is the 


Nx N matrix 


Oo 2. OO seu O =2 
=—— © FT O scx 0) 8 
0-=1 0 1 «oO 0 
a 1 ae a, tee 
mS ane — se ie 
0 0 0 90 Oo ol 
A Oe oe -1 0 é 
* 
D is anti-symmetric (so D = -D and, hence, D is 


normal) so it has eigenvalues that are either 0 or pure 
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i sin(2rkAx) /Ax for k = 0,l,..-.,N-l. Thus, the norin 


of exp(Bt) satisfies 


llexp(Bt) |] = max |lexp(iA sin(2ukA4x)t/dx|| = 1, 
O<k<N 


where we use the fact that A is symmetric so it has real 


eigenvalues, 


Kreiss Matrix Theorem 

If the approximate evolution operators Ly are not normal, 
conditions guaranteeing stability are much harder to obtain. 
The von Neumann condition (4.7) is still necessary for 


stability (why?), but it is not sufficient to ensure stability. 
One important case in which stability conditions can be found 


is for the problem studied in Example 4.1 with A no longer 
symmetric. The appropriate generalization is to assume that the 
approximation Ly has the form Ly =A ® dy, where A is a fixed 


mx m matrix (possibly not normal) and D is an N-dimensional 


N 
normal matrix. It is easy to show that 


l|exp (Lt) || = max Jexp (AVAt) || (4.11) 
N 


where AN is any of the eigenvalues of Dy - A stability 
condition for (4.11) will be obtained below. To do this, we 
generalize (4.11) and seek conditions for the Stability of a family 
of mx m matrices A(w) , where w is an arbitrary parameter. 


That is, we seek conditions such that 
max |jexp[A(w)t] || < K(t) , 
WwW 
where K(t) is a finite function of t. Once these general 


eeonditions are found, they can be specialized to give stability 


eemditions for families of the form exp (Lt) where Ly=A @ Ds with 
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Dy normal by simply choosing A(w) = Aw where w is any of 
the eigenvalues of any of the matrices Dy : 
The basic result on the stability of families of mxm 


matrices is the Kreiss matrix theorem (Kreiss 1962): 


For any family A(w) of mx m matrices, each of 
the following statements implies the next: 


(i) There exist symmetric matrices H(w) satisfying 
H(w)A(w) + A*(w)E(w) < O and 
I < H(w) , ||H(w)|| < C for some constant C. 
(ii) llexp([A(w)t]|| < C for all t20. 
(iii) (Re d) |] CAT=A(w) 72 |] x. for some constant C' 


and all A satisfying Re i’ > 0 


(iv) There exist matrices H(w) satisfying (i) with 
||H(w) |} < K(m)c'| where C' is the constant 
appearing in (iii) and K(m) depends only on 
m and not only the family A(w) . 


Observe that for a family of matrices A(w) to satisfy 
the conditions of this theorem it is necessary that all the 
eigenvalues of all the matrices have non-positive real parts. 
Otherwise there would be some w and some eigenvector a satis- 
fying \|exp [A (w) t) ul| + © as t +o violating (ii). 

The most important relation implied by this theorem is the 
implication that (iii) implies (ii) with C = K(m)Cc' That is, 
for any mx m matrix A all of whose eigenvalues have nonposi- 


tive real parts 


lJexp(At) || < K'(m) max (Re A) |] (at-a) 7+ || (4.12) 
Re \ > 0 
where K'(m) is a finite function of m. 
An elementary proof of (4.12) has recently been given 


by Laptev (1975) and improved by C. McCarthy (private communica- 


tion to G. Strang, 1975). Lapfev observes that if v>0O, then 
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t as may be proved by shifting contours in the complex plane. 

Since each entry of tovauend = is a rational function in uy 

of degree at most m, the derivatives of the real and imaginary 
parts of each entry can change sign at most 4m times when u 


increases from -» to ©. On any p-interval, say a<u<b, 


| where the real and imaginary parts of an entry in (v+ipea) 72 


are monotonic, the second mean-value theorem implies 


t 


b . 
| cos ut f(y) du = #(ay [Sinise = 


< A max|f(u)| , 
u 


for some c satisfying a<c<b where f(u) is the real or 


imaginary part of an entry in the matrix (veipea) “4, If we apply 


this kind of inequality to the right side of (4.11), it follows 
that for all i,j 


1 


co 
ipt ce SO 64m 
fon (v+ip-A) ij du t mae (v+ip-A) ij . (4.14) 


If it is true that the matrix norm has the property that 
[Bi | £ C45 for all i,j implies ||B|| < ||[c|| , then (4.14) 


implies 


64m 1 


= t 


(v+ip-A)~ (4.15) 


max 
u 


Choosing v= 1/t in (4.13-15) gives (4.12) with K'(m) = 64m. 


3 i co e vt 00 ‘ ea 
| gee as | e*rr-ay 7? aa = & | et iyeauea) ay > (4.23) 


- etnias) + £(b) feinive) - sin(ct) 


There are three important matrix norms in which 
[B,5| < C5 for all i,j implies /|[B|| < |/c||, namely 
the matrix norms induced by the Ls Lo. and L, vector 


norms. This is shown using the relations 


m 
Bil) ™ Reng oe |B, | r 
m =m 
Bil, = sup 3. B5.94 
eget det ga MOY 2 
lvl) =2 
m 
Helle = mee E ISl , 
which hold for all matrices B. In other norms. [Bi | < Ci; 
may not imply ||B|| < ||c|| but the equivalence of all matrix norms 


implies || B || < F(m) ||C|| for some finite function of the 
dimension m. Thus, (4.12) is obtained with K'(m) = 64mF(m). 

The functions K(m) appearing in statement (iv) of the 
Kreiss theorem and K'(m) appearing in (4.12) need not be equal. 
It follows from the Kreiss theorem that K'(m) < K(m) . Kreiss 
showed only that K(m) = 0(m") as m+; this is much too 
conservative. Miller & Strang (1965) showed that K(m) = 0(c™) 
as m~+o for some constant C>l. 

In the case of a normal family of matrices A(w) the con- 
ditions of the Kreiss matrix theorem are trivially satisfied: 


if the eigenvalues of A(w) have negative real parts then 


\lexp[A(w)t]|| < 1 for all t >0O and w. 


Non-Normal Approximations 


The Kreiss matrix theorem applies to approximations of the 
form Ly = A@Dy, where A is a fixed dimensional non-normal 
matrix and Dy is an N-dimensional normal matrix. This type 
of operator Ly is commonly encountered in the solution 
of initial-value problems with periodic boundary conditions. 

On the other hand, non-periodic boundary conditions usually lead 
to problems in which the non-normality affects the N-dependent 
operator Dy. When finite-difference methods are used for such 


problems, the deviation of Dy from a normal operator is frequently 


‘small’ 


“Example 4.2: Non-normality of a difference approximation to a mixed 


‘initial-boundary value problem 


A difference approximation to the mixed initial-boundary 


value problem 


u u 
CS + oe = £(x,t) (O<x<l, t>0) 


u(0,t) = 0, 
u(x,0) = g(x) 
is given by 


se a ne We = £(jh,t) (1<3<N) (4.16) 


where uy (t) = u(jh,t) anc we set Ug (t) =0 and Uyy (+) = 2uy (t) 
“Uy 1 (t). The latter condition is an extrapolation condition 
which ensures that (4.16) is a closed system of equations. This 
approximation has the matrix representation 


ase Fr UL ae i a oi i ti noms i ies " 4. 


We ee A og ee ge og a 
at oe Se OO Gk ek ee ee ae 
0 -1 0 1 . . ° . ° 0 0 0 
1 . e . e . e ° 
™N y 2h e ° e . e 
q- (OF <0) ioe -1 O 
OM 2 O20" SOS Ee es 0 -2 ' 


The departure of Ly from a normal matrix is a matrix of 
rank 1 in the lower right-hand corner. For problems of this 
kind, extensions of von Neumann stability analysis, like that 
introduced by Godunov & Ryabenkii (1963) and extended by 
Kreiss (see Kreiss & Oliger 1973), apply. 

Unfortunately, the class of semi-discrete approximations 
investigated in this monograph include problems that cannot be 
easily analyzed either by straightforward von Neumann stability 
analysis or by the Godunov-Ryabenkii or Kreiss analysis. 

In contrast to the classical problems of the numerical analysis 
of difference methods for initial-value problems, spectral 


approximations Ly are frequently not even approximately normal. 
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5. Algebraic Stability 


In this section, we develop a theory of stability and 
convergence which generalizes the classical theory discussed 
in Sec. 4. As will be shown by examples in Sects. 6-8, this 
generalized stability theory is well suited to study the con- 
vergence of spectral methods. 


A spectral approximation 


dun 
at = LiYy + fy (5.1) 


‘ to the initial-value problem Uy = Lu + f is called 


algebraically stable as N 7? @ if 


Ny os ntw8*x(t) (5.2) 


for all sufficiently large N, where r, 8s, and K(t) 
are finite for O<t<T. 

It may at first seem that the Lax-Richtmyer theorem shows 
that algebraically stable approximations cannot be convergent 
unless (5.2) holds with E-S 0) “Ee S70 % In fact, if we 
demand that the approximations converge for all u(0) and 
f(t) ain the Hilbert space 4, this conclusion is correct. 
However, it is possible for approximations that satisfy (5.2) 
with r>0 or s > 0 to converge on a dense subset of 
the Hilbert space in which the only functions for which con- 
| vergence is not obtained are highly pathological. In fact, if 


p= r+ sT > 0 but p is smaller than the order of the 


spatial truncation error of a particular solution u(x,t) , i.e. 


NP ||Lu(t) - Lyu(t) || + 0 (N + ©) (5. 3a) 


nP |]u(o) - uy, (0) |] * 0 (N > @) (5. 3b) 


nP || £(t) - £, (>| + © (N + «) (5.3c) 


for all 0 ees oR, then (4.4) and (5.2) imply that 


{A 


|| u(t) 


Uy (€) |] * 2 (N + ©) 


for Gs tet. Thus, algebraic stability implies con- 


vergence in that subspace of \ satisfying the conditions 


(5.3). If this latter subspace is large enough, an algebraic- 


ally stable method can still be very useful although it cannot 


u(0) and 


yield convergent results for all initial conditions 


forces f(t) . Since spectral methods are normally infinite- 


order accurate, algebraic stability implies convergence for 


such spectral methods. 


In the examples of algebraic stability given in Sects. 7-9, 


we find 5S i ’ s<0, and K(t) <M. In this case, 


algebraic stability implies convergence so long as (5.3) holds 


with P< ; : Thus, the approximation need not be infinite- 


order accurate to achieve convergence. However, we develop the 


general theory of algebraic stability here in the expectation 


that it will find application to spectral methods for high-order 


equations in which p may be large. 


Rake 


Our definition of algebraic stability is very similar 
to the notion of s-stability introduced by Strang (1960). 
However, Our motivation is slightly different. Strang intro- 
duced s-stability to study the convergence of time-discretized 
initial-value problems in which the norm of the evolution 
Operator grows as a power of the time step. We shall return 
to this concept when we discuss generalized stability in Sec. 9. 
Let us give an illustration of the need for a theory of 
algebraic stability. In Sec. 8, we will discuss Chebyshev 
polynomial spectral methods to solve the one-dimens inal 
wave equation u+tu, = £ (x,t) with boundary conditions 


u(-l,t) =0O. Unfortunately this problem is not well posed 


in the Chebyshev norm 


1 2 
lull? = ux) ax - 

a A-x 
© In fact, if 

1 - dl if [x] < 

u(x,0) = 
I 0 if |x| 2e, 
} 


then the solution of Uy + uy = 0, u(-1,t) = 0 at t = 1 is given 
by 


Therefore, as e > O+ , 


I|u(x,0)|]? © 


\Ju(x,1)|| 2 ~ 2 V2e 


. st aN 
so that if L= ax! 
‘ ; 1 
L u(x, 1)II 8 “se 
in fact, PeX* |) = @ for. O 4t <2, le~* || 28 


for t-> 2 , so the one-dimensional wave equation is not 


well posed in the Chebyshev norm. 


Since the finite-dimensional approximations Ly to L 


given by Galerkin, tau, and collocation approximation (see 


Sec. 2) should converge as N+” , it follows that we may 


expect 


llexp(Lyt) || 


as N?> © in the Chebyshev norm. To estimate the rate of 


divergence of llexp (Lt) || as N+ we argue that 


Chebyshev polynomials of degree at most N can resolve dis- 


tances of at most order 1/N interior to (-1,1) so we 


€ = 1/N that 


May reasonably guess on the basis of (5.4) with 


lJexp(Lyt)|]| = als i) (N+) . (5.5) 


This result is justified by the numerical results presented 


in Table 8.3. Eq. (5.5) implies that Chebyshev-spectral approximations 
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to the one-dimensional wave equation are not stable but are 
algebraically stable with r=1/4 and s = 0 in (5.2). 
Notice that algebraic stability in one norm implies 

algebraic stability in all algebraically equivalent norms. 
Thus, algebraic stability is equivalent in all of the Ly 
norms l1<p<s<o® _ because these norms are algebraically 
equivalent in N-dimensional vector spaces (i.e., they differ 
from each other only by a fixed power of N ). To show this, 
we recall that the L_ normof a vector 4 = (ayr++- ray) 


Pp 
is defined by 


N 1l/p 
Pp 
a = q la, 5 


If q=pa with O<c<a<l, then 


i=l i=1 i= 


by Holder's inequality. Therefore, for all p>l 


’ 


i 


nP tially << tlall 


N af N a 
jPO} < = lel’ a = |jal/2 ni-a/P 
i i=l p 


| 
| 
| 
| 
iq 


ea 


p 
N N 
Bi Pp 2 Pp 
Walp = 42 lagi s (42, lal) = tall 
so that 
=-1 
N lall; < lial, < Tlaly - (5.6) 


The verification of algebraic stability for spectral 


methods leads to a general problem in matrix theory. Suppose 


that Ay (N=1,2,--+,) is a one parameter family of matrices. 
We will find conditions on the members of the family such that 
exp (Ayt) is algebraically stable. We will use only the Lo 


norm since the others are equivalent to it. 


Conditions for Algebraic Stability 


Let {Ay} be a family of Nx N-~ matrices where 
Ay = 0(N*) (N+) for some finite a. A necessary 


and sufficient condition for algebraic stability 


At 
lle Ny = o (xFus*) (N + &) 


is that there exist a family {H,} of Hermitian positive- 


definite matrices such that 


ga ii 
ie tae ae a 


-1 b . 
WH Il WHyll = O(N") (N + ©), (5.7a) 


HyAy + Ayiy < c(NHY , (5. 7b) 


c(N) < dlogN (5.7c) 


for all sufficiently large N where b and d are finite 
numbers independent of N. 


To prove sufficiency we use the Lie formula 


n 
gfe Lim Ce 2t/n) (5.8) 


noo 4 


which is valid for arbitrary matrices C and D. This 


formula is proved at the end of this section. If we define 


ae 
1 2 * 
o> ace Ag, = =F Hy Ay By | 


(5.9) 


and note that 


. 1 eo 1 
exp [At] = Hy exp iy? AYHy | ae ‘ 
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it follows from the Lie formula that 


1 1 
Ayt - n 
e = limH 2 team lish ane (5.10) 


However, it follows from (5.7b) that, since C isa 


symmetric matrix, 


Also, D is an antisymmetric matrix so that 


sas | 


l}e 1 


Therefore, (5.10) gives 


ey 


ct ae 
eo* jin, || ||H 


7 N 


proving algebraic stability. 
In order to prove that the conditions (5.7) are also 


necessary for algebraic stability we define 


By = Ay - (r+1) log(N)I. 


Therefore, 


BL.t s 
lle ™ | - o(%) am ap 
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By Liapounov's theorem (Barnett & Storey 1974) there exists 


a Hermitian positive-definite matrix Hy such that 


* 
HyBy + By H 


N 


* 
HyAy + An Hy = -I + 2(r+1) log N Hy < c(N)Hy , 


where c(N) = 2(r+l) log N . In order to complete the 


proof of (5.7) we need to estimate the norms of Hy and 


ies It can be easily verified that an explicit formula 


Deny Veni ceaty gee ee ee ee eee i di he iio . a Cw a 
peeaeens SUA ee gts ge TY Ne = ’ 7 " ' " ’ es . a ian . = 4 


SOW ieee. 


so that 


JA 


-1 
2 |[Byll [tg I 


or 


2 {|B = O(N) (N + @) (5.12) 


nil 


This completes the proof of the necessity of (5.7). 


The condition for algebraic stability given in (5.7) 
implies that for every algebraically stable problem, there is 
a new norm induced by the Liapounov matrices Hy which is 
algebraically equivalent to the original norm and in which 


the problem is stable in the classical sense. 


The above result gives a method for checking numerically 
the algebraic stability of a family {Ay} of matrices satis- 


fying ||Ay|| = O(N*) as N+: 


(i) We check that the real parts of the eigenvalues 
of Ay are bounded from above by s log N ; 
otherwise, the family of matrices Ay are alge-~ 


braically unstable. 


(ii) We introduce By, = A, - (st1)log(N)I and 
compute the Liapounov matrix Hy such that 
* : 
HB, + By Hy = cl There are several numeri 


cally efficient techniques to compute Hy 


(Bartels & Stewart 1972). 


(iid) To verify algebraic stability the condition number 
of Ry must be bounded by n> for some finite b 
as N+, Noting (5.12), it is only necessary 
to verify that the eigenvalues of Hy are bounded 


from above by some finite power of N as Nv. 


This procedure is applied in Sects. 7-8 to verify algebraic 
stability of model problems, Since (5.7) gives a necessary 
and sufficient condition for algebraic stability, if these 


conditions do not hold the family of matrices A is alge- 


N 
braically unstable. 
Proof of the Lie Formula 

To prove the Lie formula (5.8) for finite dimensional 


matrices, we use the identity 


co}" [@oP [co] 
ectd we e™e” = 1e is * e"e” 
ni (se op ¢e\foge 
= i" e n n z ee" ee" 

k=0 

n 
ean £2 n-1 sak ctD CD 
lle ee" \| < ) elle DIS lle n a ee” || 
t k=0 

n-1l-k 
«(eel ' Hott) n 


cD oc 
nile ™ = e"e™} exp{(|/c]]+] |p] ]) (Q-1/n)). 


jA 
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On the other hand, 


C+D CD 
ile ate ene" || - lcp=pe}, + o( 4) (n + ©) 
2n n 
so that 
ee 
i = (.*." ee * (n + o) 


for any K > $llep-pe || + proving (5.8). 
Eq. (5.8) is also true for certain infinite dimensional 
matrices (operators). This deep result known as the Trotter 


product formula is very useful in the modern theory of 


partial differential equations. 


- 


6. Spectral Methods Using Fourier Series 

Fourier series are appropriate to solve problems 
with periodic boundary conditions. With periodic boundary 
conditions, a stable spectral method based on Fourier series 
is usually accurate and efficient. On the other hand, when 
Fourier series are used to solve non-periodic problems 
(including problems having periodic initial conditions 
but whose evolution operators violate periodicity) , 
stability is not enough to ensure convergence to the true 
solution of the problem. An example of the latter effect 
was given in Example 1.3. In this section, we investigate 
the stability and convergence of spectral methods based on 


Fourier series. 


Example 6.1: Constant-coefficient hyperbolic equation with 
periodic boundary sandr tions 


Consider the one dimensional wave equation 


uy, +u. = 0 (O< x< 1), (6.1) 


u(x,0) = £(x) 


with periodic boundary conditions 


u(0,t) = u(1,t) . 


Since collocation, Galerkin and tau methods are identical in 


the absence of essential boundary conditions (see Sec. 2), 


let us analyze the Fourier-collocation or pseudospectral 


We introduce the collocation points 


method. 
x, * n/2Nn 
a =(Ugee++eYgney?) where u, = u(x) . The collocation 


and the vector notation 


(n = 0,...,2N-1) 


equations that approximate (6.1) can be written as 


where C and D are 2N x 2N matrices whose entries are 


~ 1 - - 
kk aa exp (-2ri(k N)X,1, (6. 3a) 


Die = -2vi k' Sie F (6. 3b) 


where k' = k-N (1 < k < 2N-l) and k' = 0 if k=O. A 


simple derivation of (6.2) is obtained by observing that 


cu gives the Fourier coefficients of the collocation projection 


Pu of u(x) . Thus, pcu are the Fourier coefficients of 


- ve Pu and, finally, on} peu gives the collocation projection 


of - 7 Pu which is - P A. Pu. The matrix cC is a unitary 


1 


matrix so C* =C-, and the matrix D is skew-Hermitian so 


pD* = - D. Therefore, c7loe is skew-Hermitian so that 


llexptc7+p c)t]] =1 . (6.4) 


This proves that the Fourier-collocation method is stable for 
(6.1). The results of this example can be generalized to a 


general system of constant coefficient hyperbolic equations. 


Example 6.2: Variable-coefficient hyperbolic 


uation with 
eriodic undary con ons 


Consider the system of equations 


u, + A(x) uy = 0 0 


la 
rd 
1A 
~ 


with periodic boundary conditions u(0,t) = u(l,t) and periodic 


inhomogenity: A(x) = A(x+l) for all x. Here u(x) is 


a vector of m components and A(x) is an mxm matrix. 


If we assume that A(x) 


is a symmetric matrix and that 


for some finite a , then the Fourier-Galerkin method is 
stable. To show this, we denote by Uy the N-term Fourier- 


Galerkin approximation of wu. Then using integration by parts 


we obtain 


1 1 : 
d J * u dx = u®(AtA). udx <2a uu. dx 
at Uy Un N or = cr C* 
0 0 0 
Therefore, 
1 1 
x t = 
4 wy (x,t) uy (x, trax < & J un (x, 0) uy (x, 0) ax 
0 0 


which proves stability. 
Condition (6.5) is not sufficient to ensure stability for 


the collocation method. Consider the scalar equation (m = 1) 


U, = r(x) uy Gt eS 4 


(6.6) 
u(0,t) = u(l,t) . 


If we impose the additional restriction that r(x) is non-zero 
within 0< x < 1, then we can prove that the collocation 
is stable. We must show that exp (RC*DCt) is stable where C and 


D are given by (6.3) and R is the matrix with entries 


Ri = r(x;) 845 ‘ 


The matrix r-? can be identified as the Liapounov matrix Hy 


invoked in (5.7) and, therefore, the method is stable: 
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Rt crc*pc) + (c*p*cr*) R? = 0." 


In fact, following the proof of the main result in Sec. 5, 


Jexp(RC*pct) || * < WRIWR AH < max |r(x)|/min |r(x)|, 
O<x<l O<x<1l 


proving stability for N+o. 


If r(x) has a zero within 0<x<l, collocation with Fourier 
series may lead to instability. For example, if N = 2, the 
eigenvalues of RC’ DC are 0,0, tV"(rp+r) (ry +r3) where r; = r(x,), 
so there are growing modes if (rotrs) (r +r) <Q. In some cases, 
’ these modes may have large growth rates. One way 
to limit the growth rate of these modes 4s te rewrite 
(6.6) as 


ut 5 (r(x)u), + ; r(x)u, - ; ap Ux) u= 0 (6.7) 


Now Fourier-collocation gives the matrix equation 


> ® x + 
Uy + (&C DCR + RC DC =- Q)u = 0 
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where OQ, = - : r' (xX) ) 80 - The first two matrices on the 
right side add up to a skew-Hermitian matrix. Also, if (6.5) 
holds for r(x) then Q <* 5 aI, Therefore, we obtain the 


inequality 


Thus, we see it is possible to bound a priori the growth of modes in the 


Fourier-collocation method for variable coefficient problems with 
periodic boundary conditions. 

On the other hand, for problems with non-periodic boundary 
conditions, Fourier-spectral methods can produce wrong solutions 
even when they are stable. This is illustrated by Example 1.3 


which we now study more carefully. 


Example 6.3: 


Hyperbolic equation with non-periodic boundar 
conditions 


Consider the problem (1.7): 


mus % dul, t) ane & G<e<e te > 6 
u(0,t) = 0 (t > 0) (6.8) 
u(x,0) = 0 (O < x < mn) 
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The solution is 


u(x,t) = xt. 


If we attempt to solve (6.8) by Fourier sine series using the 


Galerkin procedure we obtain 


N 
=] a, sin nx (6.9) 
n=] 
da N 
eee nm Pe fey 4 
Oe eae et tae eS pea) 
m+n odd 


where = 0 if n is even and = 1 if n is odd. 
It is easy to verify that the above approximation is stable. 


If we write (6.10) in the form 


ea > 
ae =~ Ay até 


> 
where a = (Ayre+-r ay), f = (£).---+fy)- -. = 2 Heise 2te /n), 
then 
* 
Ay + Ay =#Q . 
Thus, | |exp (At) | | =i fer all WN and t. 


In Pigs. 6.1-6.4 we plot the solution of (6.9-10) at 
t= 1 for N = 25, 50, 75, 100 . It is apparent that Uy (x, 1) 
does not converge to the exact solution xt at t=1 as 


N+ © . Instead, u for N even appears to be converging as 


N 


at aia aa 


Fig. 6.1. A plot of uy(x,t) vs x for N=25 and t=] where uy(x,t) is 
determined by numerical integration of (6.9-10) with negligible 
time-ditferencing errors. A plot of the exact solution 

xt at t=l to (6.8) is also given. Observe the apparent 

divergence of uy(x,t) from the exact solution for 0<x<t 

and the enhanced Gibbs phenomenon at x=0,r. 


Fig. 6.2. Same as Fig. 6.1 except N=50, t#l. 


Fig. 6.3. Same as Fig. 6.1. except N=75, t=#l. 


Fig. 6.4. Same as Fig. 6.1. except N=100, t=l. | 


Vigo (2) 


-5 { 
| 
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N+ o to the function 


xt ae 
a (6.11) 
* wu (x-t) +xt st, 


for t< T, while wu. for N odd appears to converge to 


the function 


xt x > 
(6.12) 


t 
u 
odd wT (t~-x) +xt x < € 


for t <n . The results plotted in Fig. 6.5 for Uy 09 (Xr t= 2) 
are also consistent with convergence to the wrong solution 
(6.11). Notice that the approximations Uy (x, t) plotted in 
Figs. 6.1-5 all exhibit a large region of nonuniform con- 
vergence near x = 0 and x=T and that the errors in 
the interior of the interval 0< x <n decrease with N 
roughly like 1/ nN ‘ 

The origin of the divergence of (6.9-10) from the exact 


solution to (6.8) is not instability; rather, the divergence is 


due to inconsistency. Since llexp(A,t) II = 1, the method is 
stable. To show that it is not consistent we estimate the 


truncation error in the Ly norm, 
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Fig. 6.5. Same as Fig. 6.1 except ie t=2. 0 
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é, = |{tu- Lyell. 


for u = xt where Ly = PyLPy and Py is the Galerkin 
projection operator and L = - 3/3x . This error can be 
bounded from below by 


ey = ||Lu-Pybu+PyLu - PLLP.ul| 


N 


> |P yb (I-P,.) ull - |\(Z-P,) Lu|| ‘ 


However, || (I-P,) Lul| +0 (like 1/YN) as N+ because 
this norm is just the error in the Fourier sine series ex- 
pansion of Lu = - * xt = t . Therefore, if we can show 
that ||P (I-P,) u|| does not approach zero as N +o 


then (6.9-10) is not consistent. 


To estimate = ||PyL(I-P,) ul| we proceed as follows. 


Since 


(I-Py)u = ») ay (t) sin nx 
n=N+1 


we obtain 


N 
Pyyb (I-P,)u = a tie sin nx 


4 5 nm 
b(t) s-= (t) . 
n' " shin ae a =m 

mtn oda ™ ~™ 


Therefore, since the Fourier coefficients of u are given by 
n+l 


a, (t) = 2(-L) 


t/n , 


N 


2 
}PyL(I-Py)ul] = : fh, 
n=l 


2 


for suitable constants C and C)- This completes the proof 


that ||Lu - LyYl | does not approach zero as N > », 


Blair Swartz (private communication, 1976) traces the 


(6.9-10) to the incompleteness of the set 


inconsistency of 


of functions {L(sinnx) = =- ncosn x, n#l,2,...}. This set of 


functions is made complete by augmenting the set by the function l. 


Whereas u_ may be well approximated by a function Uy of 


the form (6.9), Lu may not be well approximated by the 


function Lu,. In fact, if | {Lu - Luy||+0 as No, 


then 


(Lu - Lu,) adx + 0 (Nem), 


T 1 
J Lu, = - J y na, cosnx dx = 0, 


Lu may be well approximated by Luy only if 


T 
0 = i Lu dx = u(0) - u(m), 


which is generally not true. 


As shown in Figs. 6.1-5, Wy (x,t) does converge to xt 


as N+, The analysis given above provides no clue to the 


fascinating way in which the method achieves this divergence. 


There is no indication of the 'error' wave (-2) Na (xt) that 


appears in (6.11-12) and propagates with speed 1 across 


Oo<x<m. It seems that the complete mathematical analysis 


of the divergence of (6.9-10) is difficult and we do not now 
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a 


have a justifiable argument to demonstrate convergence of wy 


to and Uoaa given by (6.11-12) as N+ @ through 


Yeven 
even and odd values, respectively. 


In the next example we will show that it is not simply 
the presence of boundary conditions but rather the non-periodic 
nature of the problem that causes the divergence of the 


Fourier-spectral methods. 


Example 6.4 Non-periodic boundary-free problem 


Consider the problem 
io + (x-5) au = 0 (8 <2 < #) 


(6.13) 
u(x,0) = £(x) 


The problem is well posed without specifying any boundary con- 


dition. However, since the solution is given by 


u(x,t) = £(5 + (xpe ) (6.14) 


it is clear that the solution is not periodic in x. Since 
r(x) =x - 5 has a bounded derivative, it follows from Example 
6.2 that Fourier-Galerkin approximation to (6.13) is stable. 
Nevertheless it is not convergent as shown by the results 


plotted in Figs. 6.6-8 for f(x) = sinx and N= 5, 10, and 


20 retained terms in the Fourier sine series. 


1 


= 


Ba i oh ha 


Pa 


u(x, .5) 


yin 
"fe 


Fig. 6.6. A plot of U(x,t) and YU. (x,t) vs x for N=5, 
t=0.5. Here u(x,t) is the exact s¥lution of (6.13) and uy(x,t) 
is the Galerkin approximation to this solution using an N term 
Fourier sine series expansion. Observe the apparent divergence 
of Uy (2, &) from u(x,t). 
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Fig. 6.7. 


Same as Fig. 6.6. except N=10, t=.5. 
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t=.5, 


except N=20, 


6.6. 


Same as Fig. 


Wags is. 


Fig. 6.8. 
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P al § ractions for Non-Periodic Problems 


There is a method that can be used to ensure that Fourier 
series yield convergent results for non-periodic problems. 
The idea is to express the solution as the sum of a low-order 
polynomial and a Fourier series; the polynomial is chosen so 
that the Fourier series converges rapidly as suggested originally 
by Lanczos (1956,1966) . The method has been used by Orszag 
(1971¢) and Wengle & Seinfeld (1977) to solve problems with 
non-periodic boundary conditions. We illustrate it here 


for the problem discussed in Example 6.4. 


Example 6.5 Polynomial subtractions applied to Fourier series 


The Fourier sine series expansion of the exact solution 
u(x,t) to (6.13) converges slowly because, in general, 
u(0,t) #0 and u(m,t) #0. This slow convergence of 
the Fourier series of the exact solution implies that Galerkin 
approximation is inconsistent, as shown using the methods 
of Example 6.3. In order to avoid slow convergence or even 
divergence, we proceed as follows. 

We seek the solution to (6,13) as the sum of a linear 


polynomial and a Fourier series: 


u(x,t) = b(t)x + c(t) (m=x) + J a(t) sin nx (6.15) 
n=1 


where b(t) and c(t) are chosen to ensure that a, (t) > 0 


rapidly as n+ © - Substituting (6.15) into (6.13) gives 
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b'(t)x + c'(t) (m=x) 


y 2nm 

m=] me mt 
n¢meven 7 ~™ 

n¥m 


a(t) = x + ‘ a, (6.17) 


are the Fourier sine coefficients of § ( 5 - x) * x an sinnx. 


, | If we knew u(0,t) and u(",t) we could set b(t)=u(",t)/7 and 


e(t)=u(0,t)/s; with this choice, the Fourier sine series in 


(6.15) does not exhibit the Gibbs phenomenon and a(t) =0(1/n?) 


as n*~, However, the boundary conditions on wu are not known 


Therefore, 


as part of the specifications of the problem (6.13). 


we must solve for b(t) and c(t) directly from the differential 


equation. 


Equating coefficients of sinnx in (6.16) gives 


= Lie wor 2 = n+1 =n 2 ‘i - 
= [c'-b'+c bla ( 1} + [b-c 2e'le e+ a, (aalecsx) (6.28) 


where “, * l1 if mn is odd, 0 if n is even; here we use the 


Fourier sine series expansion of 12 and xX: 


Pins Zarge) sap 


Hi ate LIA Dial say oe, 


oo 
ee 4 N sinnx 
_ n 
R=saa 
Ds _7)nt+l1 sin nx 
x= 2 J (-1) —— 


n=l 


Also, if b(t) and c(t) are chosen so that a =0(1/n>) 


as neo, then the Fourier series yay sinnx may be differentiated 


termwise so 


Y a, sinnx= (F-e } a, sin nx = (5- x) } na_cosnx. 
n=1 n=l) " 


n=1 


Therefore, 


hed ~~ 
F T 
lim a sinnx =95 na 
ook net n 2 } n’ 


Nis 


oo A wo 
lim /) a. sinnx = - } inl)” w €.. 
xeT- n=l - n=] ss 


Using these results and setting x =7 and x= 0 in (6.16) 


gives, respectively, 


—— 7g 


ate 


Galerkin approximation reprouuces wie eyuurces 
*, * 0 for n= N+l, Nt2, ... , 


The above derivation suggests, but does not prove, that 


a(t) + 0 sufficiently rapidly as n+ that inconsistency 


problems are avoided. The exact solution of (6.13), which 


, 3 
satisfies (6.18-20) with N=, does satisfy ._™ 0(1/n°) 
as n-o , However, the Galerkin approximation with finite N 
does not yield such a rapidly converging result. In fact, 


estimates like those given in Example 6.3 show that 


1 
\|Lv - Lyv II = wars (N > ©) (6.21) 
where v_ satisfies v(0,t) = v(1,t) = 0 and L = (Fos = 


Since the Galerkin approximation (6.18) is stable (see Example 


6.6), we expect that the errors in the Galerkin approximation 
(6.18-20) are of order N72/? for fixed ¢t. 
The above prediction has been tested numerically. In 
Table 6.1 we list for various N the maximum errors in the 
approximation obtained by solving (6.18-20). A plot of the 
error Uy (x,t) - u(x,t) vs x for N= 30, 40 at t= .5 is 


given in Fig. 6.9 =~ 10. 


it aa an ns D's 


In the next example, we prove that the method of 


. polynomial subtraction used in Example 6.5 is stable. 


Example 6.6. Proof of stability for polynomial subtractions 


It is not obvious that the approximation (6.18-20) is 
stable. Fourier series approximation without polynomial subtractions 


are stable but not consistent (see Example 6.4). On the other hand, 


Table 6.1 


4.19 
2.13 (=3) 
1.13 (=3) 
8.28 (-4) 
5.76 (-4) 
4.70 (-4) 
3.64 (-4) 
3.13 


Table 6.1. Errors in the polynomial-subtracted Fourier 


series approximation Uy (x, t) given by (6.22) and 
(6.18-20) for the problem (6.13) with f(x) = sinx 
for t=.5. Cbserve that the errors appear to decrease as 


n"3/2 as N+ © in agreement with the estimate (6.21). 
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U3q (x. -5)-ul(x, 5) 


~5x10 


Pig. 6.9. A plot of the error u ,(x,t)-u(x,t) 
in the polynomial-subtracted Fourier-series 
approximation to (6.13) with f(x)zsinx, N=30 
an@ t=.5. The approximation uy(x,t) is given 
by (6.22) and satisfies (6.18-20); the exact 
solution u(x,t) is given by (6.14). 


t 


() 


5x10 


Fig. 6.10. Same as Fig. 6.9. except 


Ugg (Xe 5) - u(x,.5) 


AA Bia t) )_f\ 


AAAS \ / 


N=40, t=. 


5. 


the approximations obtained by polynomial subtractions are consistent 
as shown by (6.21), but their stability remains to be shown. 
To demonstrate stability of (6.18-20), we reformulate 
these equations in terms of Uy (x,t) defined by 
N 


Uy (x,t) = b(t)x + c(t) (w-x) + , sw tt) einms. 
n=l 7” 


In terms of Uy (x,t), (6.18) is equivalent to 


TT . 
—e * (x >) en Sinnx dx = 0 (n=1,...N), 


w fa f) 
fs ty 
0 


while (6.19-20) become, respectively, 


Multiplying (6.23) by na, suming from n= 1 ton =N, and 


noting that 


we obtain 


% 3 3 ) 
j J =H + (x- 5) =H ft dx = 0. (6.26) 
ax 


Integrating (6.26) once by parts and using (6.24~25), we obtain 


3 3 3 
C) On zat ON - 
pete p sh] pax so. 


omar 
x 


Therefore, 


3 n fan) 2 n fou 2 1 o y (%y 2 
& GS ax-- 2] (GH fo ptegh Ee ax 


Integrating the second integral on the right once by parts gives 


n [9 2 n 3 2 « 3 2 du, \ 2 
& £8) ox-- | (32) sali (sy | « (32) 


X=T x=0 


so that 


Thus, we obtain the stability estimate 


du, (x,t) 
0x 


Nl Mai i 


The bound (6.27) shows the stability of (6.18-20). 


Examples 6.5-6 suggest that by subtracting polynomials of 


higher and higher degree from u(x,t), the residual Fourier 


series can be made to converge faster and faster. Subtracting 
a linear polynomial as in (6.15) gives Fourier approximations 
with errors of order n79/2 as Neo; subtracting a quadratic 
polynomial gives Fourier approximations with errors of order 
n’’/2, and so on. In the limit we disperse entirely with 


Fourier series and obtain a rapidly converging polynomial 


approximation. The convergence theory of these polynomial 


spectral approximations is discussed in the next two sections. 
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7. Applications of Algebraic-Stability Analysis 


The main result of Sec. 5 does not provide us with a 
systematic way of constructing the family Hy of Liapounov 
matrices necessary to prove algebraic stability. In general, 
these matrices are difficult to find. However, there are 
several problems for which they can be found directly from 
the differential equation, 


It is very easy to construct Liapounov matrices for Galer- 


kin approximations to 


Ye ty 


where L is a semi-bounded operator on the Hilbert space x. 


We say that L is semi-~bounded if 


L+bL < af (7.1) 


for some constant a , Where L* is the adjoint of L defined 
with respect to the Hilbert space inner product (,). If L 


is semi-bounded 
Aiu,u) < a(u,u) (7.2) 
at , Pa , e >. 
so 


(u(t) ,u(t)) £ e**(u(o) ,u(0) 
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and the ‘energy’ (u(t),u(t)) grows at most exponentially with 
t . 
If an energy estimate of the form (7.2) exists, then Galerkin 


approximation based on the Hilbert space inner product (°,°) is 


stable (and, hence, algebraically stable). The Liapounov 
matrix H, may be chosen to be the N xN _ identity matrix 


In - In fact, it follows from the Galerkin equations (2.6-7) that, 
if f = 0, then 


ar (uy Uy) = (Uys (L4L") uy) S$ aluy,uy) 
re 
i 
Thus, 


(uy (t) pay (t))¢ e%* (uy (0) ,uy (0)) 


Since Uy (t) = exp (Lyt) uy (0) for all Uy (0), it follows that 
llexp (Lt) || < exp(#at) so stability is proved. phe reader is 
reminded that with stability established, the theory of Sections 


4 and 5 proves convergence for consistent schemes. 


| Example 7.1: Semi-bounded Galerkin approximations 
The above construction establishes stability and thus con- 
vergence for a wide variety of Galerkin approximations. Among 


these stable Galerkin approximations are: 


(i) Solution of any problem % = Lu that is semi-bounced 


(-1,1) by means of Legendre series. For example, 


} in Ly 


utu = f(x,t) with u(-l,t) = 0 is stable (and convergent) 


when solved by Legendre-Galerkin approximation. For our argument 


to be complete it is necessary to verify that the Legendre- 
Galerkin approximation to this problem is consistent. This is 
done as follows. 


We write 


|| Lu-P LP ul| < || (I-Py) Lul| + [|P i (I-P,) ull . 


N 


The first term on the right goes to zero as N*o~ ata rate 
governed solely by the smoothness of Lu ; it measures the 
error in the N term  lLegendre-Galerkin expansion of Lu. 


The second term is estimated as fullows. Set 


L(I-Py)u = as ano, (*) 


where {¢} are normalized Legendre polynomials. If L is 
a finite-order differential operator so L* is also a finite- 


order differential operator (for example, L*=3/3x if L=~-3/9x), 


then 
an = (oo. L(I-P,) ¥) 
= (Lo. (I-P,) ¥) ° 
Thus, 
x 
lanl < He ott [l(r-py el] 


O(n" /n®) — (nv@ 7 Nem ), 
where A depends only on L (A = 3/2 if L = -3/ax and on is a 
normalized Legendre polynomial) and B depends only on the 


smoothness of u (B is arbitrary if u is infinitely differentiable). 


Thus, 


||Pyb(1-P,) ul > 0 


faster than any power of 1/N if u and all its derivatives 
are smooth. This proves consistency. This kind of proof 
extends to a wide variety of the examples to be discussed in 


Sects. 7 and 8, but will not be repeated. 


(ii) Solution of u, = xu, with the boundary conditions 


u(tl,t) = 0 is a well posed problem in the Chebyshev inner 


product : 
(uv) = f BGdvOd ay , 
-1 (1-x4)# 


In fact, if L =x 9/ox , and u is differentiable and 


satisfies u(tl) = 0 then, by integration by parts, 


1A 
o 


1 or 
{ (-x?)7 Fue ax 


(u,Lu) = f? x(1-x?)7#u ujax = - 
ee, | “I 


Thus, Galerkin approximation to the problem is stable using 
Chebyshev polynomials. 


(iii) Solution of u, tu, =0 (0 < x < =) with 


u(0,t) = 0 is a well posed problem in the Laguerre inner 


product 
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(u,v) = f ulx)v(x)e *ax . 
0 


In fact, if u(0,t) = 0 then, by integrating by parts, 


@ eo 
ae | uu e “dx = - $ emer it at e*u*ax < 0. 
0 0 


Similarly, the problem u, = u,, (0 < x < ~) with u(0,t) = 0 


is also stable in the Laguerre norm. 


(iv) Solution of ur = ~xu, (-~ <x < ~) is well 


posed in the Hermite inner product 


o 2 
(uyv) = f ulx)v(xde™ ax . 


—o 


In fact, 


3 ro 
sete) = -2 Jfxe 


oO 


2 
* uu, ax 


so that integration by parts gives 


x 


© 2 
9 (u,u) = f u2e* (1-2x)ax < (u,u) 
st = 


—a@ 
where we assume that u << vx exp 5 x?) as |x| + 
(v) The heat equation % = a. with u(+l,t) = 0 is semi- 


bounded in the Chebyshev norm. In fact, if u is differentiable 


for |x| <1 then 
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1 =- - 1 -i 
| (1-x?) 2 uu, dx = (1-x?) uu | - | fatiex®) 2) u. dx 
ae -} x x 


eae 
* fee 


XxX 
* ek 


The first term vanishes because u is a polynomial in x and 


a teat 32 


{ therefore u(tl) = 0 implies 


The integral term on the right is 


1 2 2 
- | (u(1l-x”™) ) we ax 
at x x 


1 
2 2 2 
= - [are ) ) , (u(lox ) ) , (1-x ) ax 


Nie 


er -i 
1 1 9 2 2 
+s bs 5x | (Ulex) x(1-x") dx 
1 1 
1 ne . aoe 
= - | (u(l-x") ) (l-x") ax 
-l ' 
2 «ik 2 2 
# 2a eta-0") -% a? (1-4) ax (7.3) 
2 Zz 
-l -1 
< 0 
and therefore 
1 2 
= [ a axe <¢ ¥ . 


-1 Vi-x? 


In the next two examples we generalize the proofs of 
a. stability and convergence for Galerkin approximations given in Example 


7.1 to show the stability and convergence of tau approximations. 


Example 7.2: Semi-bounded tau approximations 


(z) Cohsider the equation 


du du 
gt ox 


with 


u(tl,t) = 0 


It was shown in Example 7.1(ii) that if L= x3/ax, then ' 


b+ Lf <0 
in the Chebyshev inner product. If we seek the solution as the 


truncated Chebyshev series 


{ N 
On * i s Ty 
i n=0 
by the tau method, then uy, satisfies exactly the equation 


9 5 

q ed = “Ny = Tyg) Ty 00) + Tynz (t) Ty-} (*) (7.4) 
q ot. ox 

i 

i} 

q Equating coefficients of x ana xN-l on both sides of 


i | 

q (7.4), we obtain | 
| 
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ae 


gn-1yn ma n2n73,n-2 


since * = + +*** . Therefore, 


| ay 
| (0-3) = RGR Ugg) Fy Mag ay 
1 


conditions u(+l,t) = 0. 


e - — 
+ fayey7 (N-L) ayy) ay-y (7.5) 
| so that 
{ 
<s uu oe ed iy on . a <4 
2 5E fiuy ay) ay inl = ((L+E)uy uy) Nay - (N-l)ay_, = 0 
Since 
N 
2 
(UyeUy) = YL ay 
n=0 
i 
| the above inequality is equivalent to 
N-2 
| 3 ) 2< 
a. = 0 (7.6) 
vw neg COS 
This proves stability: an and Aay., are bounded because they 
are determined in terms of Agr Ayre++rAy_ by the boundary 


For this example, we can prove stability directly from the 


i matrix representation of L,. In fact, 


| 
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1 N=4 
(au We Pek + 2 RES ean ek 2 Ie 
mege Cyr aR pS K 4th, pI COSISN, O<k<N), (7.72) 
2 even 
In the tau approximation, the boundary conditions u(tl,t) = 0 require 


that the last two rows of the matrix Ly be replaced by 


k 
= (-1) (7. 7b) 
yet, k : 
(L..) dan (7.7c) 
Boy 


If the boundary conditions (7.7b,c) are not applied then 
the spectral approximation is unstable: without the boundary conditions 
N 


Ly has the eigenvalue N_~ [with the eigenvector Ay-2K7 'q) 


aye2K-1 = 0) so that 


To prove convergence when the boundary conditions (7.7b,c) are 


applied, let us first consider an oda solution in which oe 0 


if n ‘$s even. If we assume that N = 2Mtl and set 


Gy = F2K+1 (cof k <M) 


then the system reduces to 


Mend 
ot 
where 
M-j 
D5y =~ (2K41) 55,4 2 ce (2K4I) 8 549 pr 2N COSGS My O<k<M) 


If we introduce the M x M transformation matrix S defined by 


S = 6 = 
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then S(D+D*)S* is a diagonal matrix with entries 


(-4, -4, ..., “4, ~4N - 12). Thus, we obtain D+ D* < 0, 
so that a(d,d) /at<0 which proves stability. 


Example 7.3. Stability of tau _ methods applied to degree-reducing 
semi~bounded equations 


An argument similar to that given in Example 7.2 demonstrates 
stability of tau methods in terms of arbitrary orthonormal polynomial 
a = Lu where L is semi-bounded and degree 
reducing: L is said to be degree reducing if for any polynomial 


bases for equations 


Py of degree N, LPy is a polynomial of degree at most N ~- k 


where k is the number of boundary conditions that are applied. 


If L is degree reducing, equating coefficients of x 7 peverX 
in 


dU, 
bike nee key 7 2 


implies that T(t) = ay (t) for n = N-k+tl,...,N; here 


N 
Uy (Xt) a 


-£, a(t) 6, 0) 


and the orthonormal expansion polynomial $9) is assumed of degree n. 


Therefore, 


ra 
4 se (uy dy) - x ata, = ((L+L*Ju,,u,) § 0 


so that 


eae 


c 


t 
4} 
’ 


which proves stability since Ayiypyyrees Ay are determined by 


the boundary conditions in terms Of Apr Ayre+-rAniy ° 


Example 7.3: More stable tau approximations 


(i) Suppose that 


m+ aS 0 (el < x ££ 1, t ? 0) 


u(-l,t) = 0 


is solved by tau approximation using Legendre polynomials. 


The hi degree Legendre polynomial Uy satisfies 


so that 


1 
a 2 2 
a '{ No - mar ay! < 0 
which proves stability. 
(ii) Suppose that 
~— Yxx 
u(+1,t) = 0 


is solved by the tau method using Chebyshev polynomials. Since 


2 * 
L = +, is degree decreasing and L+L *0 (see txample 
ax 
7.1(v)), the method is stable. 


(iii) The solution of 
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ee 


pee DU kn 


= sint (t > 9) (7.8) 


c 

Oo 

ct 
i 


c 
* 
2 
oO 
S 
“A 
* 
A 
= 


by Laguerre polynomials is stable using the tau method since, 
by Example 7.1 (iii), L is semi-bounded. The equations of 
the Laguerre-tau approximation to (7.8) are a simple modification 
of (2.23-24). In Fig. 7.1 we compare this tau approximation 
with the exact solution of (7.8) at t = 30 for a 20-term 
Laguerre expansion. The reader should compare this approximate 
result obtained by the tau methoec with the best Laguerre approxi- 
mation to sin x plotted in Fig. 3.12. 

In the next example we discuss some ways to find non-trivial 


Liapounov matrices {Hy} when L is not semi-bounced. 


Example 7.4: Polynomial approximations to a_ variable coefficient _ 
hyperbolic equation 


Consider the initial-value problem 
= -xu ix{ <2 
u(x,0) = g(x) (7.9) 


which is well posed without requiring any boundary conditions. 


The exact solution.wto this problem is 


t 


u(x,t) = g(xe~ 


) 


-— (oe x) 9%n 


(X-0€ ) wi 


"ete *bta 

uy pezzoTd (o¢’x)N 03 UOTIeUTxOIdde eizenbey 3s0q 

ay? OF ucoTJeEWTxOidde nejz-erztEeNhbey syQ jo AZtzerttwts 
=U 

eu} ePATeSqO “UTS = Ue "3 Aq peoetdez (#7z°z) 

pue 0=3 UTA (€7°z) Aq usatb oze suotzenbse neq syQ 

pue 0Z=N 3e& pezeounz3 st (zz‘°z) uotsuredxea ezrzsnbeq7 

243 220H “OZ=N ‘O€=3 I0Z weTqoad sty OF UOTIANTOS 

3oexe ay3 Jo 30Td e pue (g°Z) weTqozd ay3 03 (3/x)Nn 


UOTIBUTxXOIGae Nej-sizenbey ey3 30 30Td WwW T°L “btaZ 


— 


so that u(x, t) approaches a constant as t?*™; 


u(x,t) ~g(0) (t+, |x|<1). 


The problem is well-posed in the sense that ||exp(Lt) | | 


is finite for finite t, where L= - x3/3x and Ill] is 
any t because the function that extremizes | |u(t)|| subject 


the usual L, norm. However, [lexp (Lt)|]= exp (5 t) for 
to {|u(0)||= 1 satisfies u(x,0)= g,{x) where 
{ 


t/2 P 
+ 2 |x\<e"* 
g, (x) = be 
. 0 |x|>e7*. 
Therefore, ||exp(Lt)|| grows exponentially as t>>. 


The operator L is semibounded in the usual Ly norm: 


1 1 2 1 

r) 2 du: 2 2 2 
ara uo «dx = - X wo dx = - u' (1) - ui (-1) + u ax 
ve Ee ax e 


1 
< J u” dx, 
a) 


so L + L* < I. Therefore, Galerkin polynomial solution of 
(7.9) is stable and convergent. The Legendre polynomial approx- 
imation Uy (x, t) satisfies 


(7.10) 


—— > 


3 


dUy dUy 
t 


* 
ow 
| 

" 

oO 


aS 


a ae 


oe ey ERNE oar a 


exactly because no boundary conditions are applied and L is 
degree preserving. Therefore, Galerkin, tau, and collocation 
approximations to (7.9) are identical and all three methods 
are stable. 

In fact, all polynomial-spectral methods applied to 
(7.9) satisfy (7.10); all polynomial methods for this problem 
give identical results and, therefore, they are all stable in 
the usual Ly norm. In terms of the natural norms for a general 
polynomial basis {fy}. i.e. that norm in which (p55) =8 55" 
the spectral approximation (7.10) is algebraically stable if 


the N x N' matrix whose elements are 


(Hy) 5% = f v5 (x) Wy, (x) ax 
has a condition number which is bounded algebraically, i.e., 

-1 B 
HE agl tf lay 71] = 0 cn8) (nee), 

As an example of the complicated behavior of spectral 
approximations for this problem in norms different from the usual 
Ly norm, let us consider the Chebyshev-L, norm. It may easily 
be shown that L + L* is not semibounded in the Chebyshev inner 


product. For example, consider the trial function 


v= To - Ton 


then 
! 


| ((L+L")v,v) = - (xv, Vv) = (v,xv,) 
ai ely 
So 


-(-anttgy-1*- --+T,1,T) - 


5 N(v,v) . 


Nevertheless, Chebyshev approximation to this problem is 


algebraically stable. This fact may be explicitly demonstrated 


by construction of a Liapounov matrix. 


' A Liapounov matrix for the Chebyshev approximation to (7.9) 


may be found by direct examination of the evolution equation for 


= . 
the vector ay = (Aagr-++ray)? 
dan N 
oc She, ~ 2 ) Pa, (n = 0,...,N). (7.11) 
p=n+2 
pPtn even 


Since a, decouples from AyreeerGy in (7-11),we can restrict 


attention to a),...,a)- Suppose we define {H,} by 


(H,) aa , (1 < 3-k < N), 


Then 


0 -1 0 -1 0 
T 
HyLy + Liaky =f -l 0 -l 0 -1 ; 
0 -1 0 -1 0 


the matrix displayed above has rank 2 and the nonzero eigenvalues 


are -([N/2], -((N+l)/2]. Therefore, by the theory of Sec. 5, 
Lyt = 
We ycVilniling i= ve, 


where ||:|| is now the Chebyshev norm. Thus, L, is 
algebraically stable in the Chebyshev norm even though Ly thy 
is unbounded in this norm. 

The qualitative behavior of | lexp(Lt) | | as a function 
of N and t is as follows. For fixed t and N+a, 
| Jexp (Lt) | | = 0 (nt/4), this result is justified heuristically 
by following the argument given in Sec. 5 that led to (5.4). 
On the other hand if t>inN, |jexp(L,t)|| = o(n?7?) as New, 
A heuristic justification of this result is as follows. Let 
u(x,0) = 1 for |x|/<e, 0 for |x|>e. Then the exact solution 
of (7.9) for t>&n l/e is u(x,t) ¥ 1 for |x]<1, so JJux,t) {]%aq 
as e*Ot for t>&n l/fe. As in Sec. 5, we conclude that 


| |exp(Lyt) | |= 0 (nt/?) for t>@nN as Neo, (Even in the usual 


L., norm, | | exp (Lt) || = 0 (nl/*) when t >&nN, which mimics the 
unbounded growth of = || exp(Lt)| | as t +m.) 
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8. Constant Coefficient Hyperbolic Equations 


In this Section, we discuss the stability of spectral methods 
for the problem 


au 4 22 = 0 isl < iy £70) 18.0) 
with the initial condition 
u(x,Q) = £(x) (|x| < 1) (8.2) 
and the boundary condition 
u(-1,t) = 0 (¢> 0} . (8.3) 


The results for this problem canibe extended to a general 


hyperbolic system of the form 


with characteristic boundary conditions, because for any hyperbolic 


system A can be diagonalized by a real similarity tranformation. 
The Operator L=- Ss is semi-bounded in the usual 


L5(-1,1) norm when operating on the subspace of 


functions v_ that 
satisfy the boundary condition v(-1,t) 


0. In fact 


1 
,(L+L* = « dv 
(v, (L+L*] v) A V 5y ax 


~v*(1) 


“A 
o 


and therefore Galerkin and tau methods are stable using 


Legendre polynomials. 


However, L is not semi-bounded in the Chebyshev norm. To 


show this, we set 


v(x) = Toy (*) - T, (x) - 2T,) (x) 


so that v(-l) = 0. In this case, using the result 


Toy 7 2NIToy-a*Tan-3t °° tT! 


we obtain 


1 
(v,[L+L*)v) = -2 f 
=e 


1 4 ad 
= -2 f (1-x*) "2 [2n(t 


x +. -+T)) -To] (Toy-T,-2To) ax 


2N-1*? 2n-3 


Ane wv). (8.4) 
The fact L + L* is not semi-bounded is consistent with the fact 
that exp(Lt) is not a bounded operator for t<2 in the Chebyshev 
norm (see Sec. 5). However, these results do not prove that . 
Chebyshev-spectral approximation to (8.1-3) is not convergent. 
In fact, we shall show that, while Chebyshev-spectral approximation 
to (8.1-3) is not stable in the Chebyshev Ly norm, it is algebraically 


stable in this norm. 


In order to investigate algebraic stability, we must wisatiel more 


carefully the behavior of the Chebyshev coefficients of the 


@pproximate solution 
Uy = LY a, (tT (xy. 


The cifferential equations for the a's are given by (2.11) 


for Galerkin approximation, (2.19) for the tau method, and 


(2.32) for the collocation method. As remarked in Sec. 2, all 


these equations may be written in the vector form 


where a= (Agray,...ay) and Ly is an (N+l) x (N+1) matrix. 
N 


umerical Evidence for Algebraic Stabilit 
——o i eegebraic Stability 


Let us first examine the behavior of L. + Ly*- In Table 8.1 


we list the largest eigenvalue of Ly + Ly* for N= 10.20)... «L106 


for the three Chebyshev methods. This table indicates that the 


largest positive eigenvalue of Ly + Ly* grows like cn? for some 
constant C. If Ly were a norma). matrix this would imply 

Lt 
that |le N || 


ee like exp (5 cn? t). However, the matrices 


Ly are not normal and therefore the large eigenvalues of Ly 


+ Ly* 
do not imply instability. 
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Table 8.1. 


The largest positive eiqenvalue \ 


max 


PIRES, — 
Collocation Tau Galerkin 

: | 68.84125 | 21.4089 72.8947 
| 287.6920. | 84.8970 296.3027 
j 30 656.6218 | 190.4908 669.6434 
| 40 | - 13.75.2228 338.1769 |  1192.9232 
| 50 | yg43.8839 | 527.9525 | 1866.1433 
60 | 2662. 4966 | 759.8167 | -2689.. 3042 
7 70 3631.0503 | 1033. 7690 | 3662. 4061 
i 80 | «749.5453 1349.8093 | —4785.4489 
| 90 6017.9812 | -1707..9375 6058. 4329 

| { aase.3se8 |= (208.894 7481. 3579 


of Ly 


for the Chebyshev-spectral solution of the one-dimensional wave 

equation (8.1-3). The Galerkin approximation to this problem is 

given by the solution to (2.11), the tau approximation is given 

by (2.19), and the collocation approximation is given by (2.32). 

ec * 0.75 for the 
for the tau 


Observe that \ \ on? as N +o where 


max 
Galerkin and collocation methods and c # 0.21 


method, 


otic ssa atirlnserfibty ate Us 


In Table 8.2, we give the norms of the matrices 
exp IL) . expiLy | for the three projection methods (Galerkin, 
collocation, and tau). Theresults indicate that | Jexp(L.) || grows 
only like ni/4 as New (as argued heuristically in Sec.5). In other 
words, L. is algebraically stable (at least for t=l). This result 
shows the extreme pessimism of the energy estimate | }exp (Ly)| | = 
0 (exp (5 cn’) ); crude energy methods may be very misleading for non- 
normal evolution operators. 

In order to understand better how the Chebyshev spectral 
methods avoid an energy ‘catastrophe' [energy growth like exp (Nt) ] 
we have solved the tau equations (2.19) numerically with a very 


"bad' initial condition: 


Uy (x,0) = (T(x) + 2 (x) + (-2)™ TQ Q01/ VT . (8.5) 


For the tau method, this initial condition satisfies 


R) 2 
yt (Uys By) =(Uy, (L th ) uy) = O(N”) (N + o ) 
t= 


0 


In Figs. 8.1-2 we plot the energy (Uy,/Uy) vs t for N= 25 
and N= 50. It is apparent that the initial slope of the energy 
growth is of order n? but that the energy does not maintain this 


rapid rate of growth. Observe that the region of rapid growth 
is closer to t= 0 for N= 50 than for N = 25, The behavior 
observed fn Figs. 8.1-2 is not inconsistent with the fact that 


uy (t = 0) is a 'bad' eigenmode of Ly + Ly*s Because Ly is 
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Table 8.2 


Galerkin 


2.0003 2.5788 


2.8119 3.1903 


3.4857 3.8328 


4.0514 4.4078 
4.5339 4.8630 
4.9855 5.2057 
5.4002 5.5262 
5.7770 5.8689 
6.1401 6.2526 


6.4831 6.6257 


Table 8.2. The largest eigenvalue Nase of exp (Ing) exp (Ly) « 
Observe that Amex behaves as oni /2 as N + © where 

c # 0.6 for all three spectral methods. The 

largest eigenvalue of exp (Ly) exp (Ly) grows only like 
ni/2 despite the existence of eigenvalues of L,, + L* 


2 N N 
growing like N° (see Table 8.1). 
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Fig. 8.1. A plot of the energy (uy), Uy) ws tt. £05 
the 'bad' initial conditions (8.5) with N = 25. 
Here uy(x,t) is the solution to the Chebyshev-tau 
equation (2.19) for Uy += 0 with N=25, 
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non-normal the 'bad' initial condition is not an eigenmode of 

Ly so that after evolution from 0 to t exp (Lt) , Uy ‘rotates' 

out of the region of bad modes of Ly + Ly rs 
The direct computation of exp [Lyt] for t=l is not enough to 

verify algebraic stability because the theory of Sec. 5 shows 

that we must study the behavior of exp [Lt] for a complete 

time interval 0<t <T . This may be done using the method 

suggested in Sec. 5 for the numerical verification of algebraic 

stability. First, in Table 8.3 we list the numerical]ly computed 

eigenvalues of L, . Observe that all the eigenvalues of 1, have 


negative real part. (This result will be shown rigorously later.) 
Therefore, [lexp (Lt) || >0 as t*= for fixed N. Thus 


the Chebyshev approximations are asymptotically stable in the 
sense that they remain bounded as t>*@ with N_ fixed, 

In Figs. 8.3-5, we plot the L, -matrix norm of exp (Lyt) 
ve t for N=5,15,25. Observe that as t+ for fixed N, 


1/2 


| |exp (Lyt) | 1) approaches zero while it grows slowly (like N°’ “) 


as N*~ for fixed t<2 (Note that growth of [exp (byt) 11, 


like yi/2 as WN? is not inconsistent with growth of 
| Jexp(Lyt) ||, like nt/4 Also observe that the norms seem 
to have a boundary layer at t=2 such that [exp (Lyt) 11, +00 


as N*o for t<2 and *0 as N*o for t>+2. This behavior 


is consistent with the unboundedness of exp(Lt) for t<2 (see (5.4)]. 


Asymptotic stability does not prove stability because Ly is 


not normal. The next step in the computational proof of stability 


is to compute numerically the Liapounov matrices Hy satisfying 
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Table 8.3 


eet ee rn eee ween ees em. 


Collocation Tau Galerkin 


-2.4532 ~2.999 -1.9306 


-2.5932 ~3.9320 -2.15 

-2.7267 ~4.5380 -2.32 

-2.849 ~4.9918 -2.4659 
-2.966 -5.3837 -2.5965 
-3.0824 -5.7266 -2.7226 
-3.1985 -6.0489 -2.8478 
-3.3162 -6. 3650 -2.9738 
-3.4365 -6.6861 | -3.1017 
-3.5597 -7.0229 -3.4335 


Table 8.3. The real part of the eigenvalue of Ly with 
least negative real part for the collocation, tau, and 
Galerkin spectral approximations to (8.1.3). Since all the 


eigenvalues of L have negative real parts, these spectral 


N 
methods are asymptotically stable as t +, 


Fig. 8.3 A plot of the L, matrix norm of exp (Lt) vs t 
where Ly is the Chebyshev-tau approximation to (8.1-3). 


Here N=5 so the Chebyshev polynomial expanison is truncated 
after T, (x). 
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|]exp(L,.t) 11, 


Fig. 


8.5. 


Same as Fig. 


8.3. 


except N=25. 
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* 
Hy” +A. HH. = = I (8.6) 


A good method to compute Hy is described by Bartels & Stewart 


(1974). In Table 8.4 we list the condition number of Hy for 


the Galerkin, collocation and tau methods. This table suggests 


that the condition number of Hy grows at most like n> as N 


for the Galerkin and collocation utioes and like n? for the 


tau method. Recalling (5.11), we obtain 


3 
[| fexpiytl| |= 0 (Ne *) . (8.7) 


1 
Nicer 


for all three methods. It should be noted that (8.7) gives only an 
upper bound for |lexp[L,t]]| . According to the theory given in 
Sec. 5, this upper bound can be sharpened by at most LILytl = 0(N7) 
(Nem), explaining the origin of the difference between the estimate 
(8.7) and the observed behavior yi/4 of the computed Lj-matrix 
norms. 
In the above discussion, we have given numerical evidence 
for algebraic stability of the Chebyshev-spectral methods for 
(8.1). We shall ‘now prove rigorously that Chebyshev-spectral 
methods for (8.1) are algebraically stable. 
_Proof of Algebraic Stability for Chebyshev-Galerkin Approximation _ 
In the Chebyshev-Galerkin approximation to (8.1), we represent 


the spectral approximation Uy by the series 


N 
aa » oo fed \ im 
un, = } a(t) [T,-(-1) "T)] (8.8) 
n=] 
+ ar : 5/2 
The condition number of Hy can grow no faster than N as 
- ie) 
Noo, To see this, wenote that (5.14) gives hte ie 0(N%) 
while (5.13) and the results that || exp(I,t)|| = O(N 4) for t<2 
and [lexp(iyt)||>0 as N°" for t>2 give |Inyll= o(n)/*) 


as N°?S 


I 
_ 


Table 8.4 


Pe PMS so SERENE TAO Ie AO FR 
N Collocation Tau a Galerkin | 


ag aren eer nes eee rn en at en ee re en ee an et ee nnn ce 


\ 
10 4.1463 x 107 ~—- 3.1090 x 107 4.6388 x 10° | 
| | | 
20 3.0332 x20" | 2.2421 2 20> | 3.9672 » 20° 
30 9.8746 x 10> : 2.7938 x 10° | 1.0464 x 104 
‘ ‘ ” { 
! { 
40 2.2940 x 104 | 4.9662 x 10° | 2.9083 x 104 | 
50 4.4220 x204 | 7.7593 x10? | 4.6138 x 104 
coder eee (SE Bie 
Table 8.4. The condition number Wig ELLER | in the 
Ly matrix norm of the Liapounov matrices Hy for the 


collocation , tau, and Galerkin spectral methods for (8.1~3). 
For the collocation and Galerkin methods, the condition 
number seems to grow like "3 as N-o , while for the tau 


method it seems to grow like n? as Neo , 


Recalling (2.34), Uy satisfies 


(~1) . (8.9) 


We can determine T(t) by equating the coefficients of x 
in (8.8): 
day (t) 
e we 
Try (t) eee tek) 


Let us now multiply both sides of (8.9) by 2(1-x) uy and integrate 


with respect to the Chebyshev weight function (i-x2) 1/2, gus, 
the left hand side of (8.9) becomes 
1 r) r) 
2] (1-2) ugl se + A (Hx 7? ax 
-1 
a; 2-4 2 : 3 -4 Py 
=a¢ [ (1-x) (1-x") ugax+ f (1-x) 7(1+x)* 4 ax 
-l =-1 
naa = ee yee: 
+ f Uy [4 (1-x) * (1+x) * + 4(1l-x) 7(1+x) 2] ax (8.10) 
-1 


«1597 
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The boundary term in the last expression vanishes because Uy 


is a polynomial satisfying Uy (1) =o. Also, 
N ny N i 
(1-x)uy = Cham) 1 ait, tt) 9) = Aqlt- ten) To) 


; n 
+ tone a. (-1) T)] 


eT 


N N . 
= T -(-1)"7,)- SP th e Hh oes 
ae yan! n (=) T9] #1 Pn ut (-1) T,) anita (-1) T,) 


(8.11) 


The first and third sums on the right in (8.11) are orthogonal 


to the right side of (8.9). The inner product of (1-x)u, 


with the second sum on the right in (8.9) gives 


da 
K "ree oe 
~#a, ge <7 ge on (8.12) 


~(-1)N 
(-1) Ty eae 


N 


Combining (8.10) and (8.12), we obtain 


v a 
sa if (1-x) (1-x?) us dx +3 & ay 


nN 
1A 


0 (8.13) 


— 


; This inequality proves that Uy is stable in the new norm defined 


: in (8.13): 


1 
Jul]? =f (rex) (dex?) 727? Jue] 2 ax (8.14) 
a 


It remains to prove that the norm defined by (8.14) is 


algebraically equivalent to the usual Chebyshev-L, norm. That 


is, we must show the existence of two functions c, (N) and co (N) 


such that for every Nth degree polynomial Uy 


2 
1 u i 
cf Me axs J —ae oxce, sf She ax (8.15) 
Ay (Vieni ~ | OF bes al fies 


where 1/c, (N) and c, (N) grow at most algebraically as Neo, 


The second inequality in (8.15) holds with ce, (N) = 2 because 


1-x<2. 


The first inequality in (8.15) is more difficult to 


establish. By the mean-value theorem, 


‘doesn 2 Roe. 
J. pie OF w. UPeE. 5 f ee Ox Lee) 
by fee Wa Joe : 


However this does not prove the required inequality because it is 


not clear that 1/ (1-Ey) is bounded algebraically as N° for all 


polynomials. 


To establish the first inequality in (8.15) we use a different 


approach. 


We substitute the Chebyshev polynomial expansion 


ae e SS 


{ et te tee ~ See + Ta * 
a a n=] 2® 
1-x 


Wey eet kee 


where a is the symmetric, positive definite, (N+l) x (N+1) 


tridiagonal matrix whose elements are 


cy if j =k 
tn “ke, if j = k-l 
4 pe (8.16) 
nN’ 3 — ey, if j = k+l 
0 otherwise, 


where cy) = 2, c, = 1 if n > 0. Yo complete the demonstration of 
the first inequality in (8.15), we must show that Hy>cy (NYT 
where ¢,(N)>0 and 1/c, (N) is bounded algebraically as N+. 
Since By is nearly a constant-diagonal tridiagonal matrix, the 
eigenvalues of Hy can be studied by standard techniques: if 


Dy = det (H\-AI). then Dy satisfies the three-term recurrence 


relation 


lle en 


Since (8.17) has constant coefficients, it is easy to solve 


exactly. From this solution, it is not hard to show that the 


smallest eiqenvalue of H 


satisfies 
2 
r (N) bi acing (N>©). 
min 8N 
Choosing c,(N) = A_,_(N) gives l/c, (N) ~ 8N7/n?  (N+>) 
- 1 min : Ve ; : 


r, This proves that the norm defined by (8.14) is algebraically 


equivalent to the Chebyshev norim and, therefore, Chebyshev-Galerkin 


approximation to (8.1) is algebraically stable. Note also that (8.13) 


shows that the Matrix Hy defined in (8.16) 


c(N) = 0. Since |/Hy]] = 0(1) ana 


satisfies (5.7b) with 


Jue | = o(N?), (5.22) 


implies that | |exp (Lt) || = 0(N) as Ne», which also follows 


directly from (8.15). 


We have not yet been able to obtain a rigorous demonstration 


that [ [exp (Lt) | | = @ in?/4) as Neo 


© 


as found numerically in 


Table 8.2. Our best result to date is | Jexp(nyt) || = O(N) as 


Noo, 


Although the problem (8.1) is not well posed in the Chebyshev norm 


(as shown in Sec. 5), it is well posed in the norm defined by (8.14). 


Using (8.1) and (8.3), we obtain 


1 
ek f.  BF ioe cen Ree 
2-1 
Thus, 
1 1/2 
d 1l=-x 2 
at I (153) ge REM 


so that |le%t||<1 in the norm (8.14). 


Proof of Algebraic Stability for Chebyshev-Tau Approximation 


The proof of algebraic stability for the tau method is similar 
to that just given for Galerkin approximation. The Chebyshev-tau 


approximation uy satisfies 


r) du 
be + Ne nt) ty) (8.18) 


Uy (-1,t) = 0 ' 


where 


(8.19) 


| 
| 


Therefore , 


j 2 N-1 

97u da 20 

9 OE cc ee A (8.20) 
(1-x) Jxot =-N at Ty wtghh n 


Moreover, comparing the coefficients of - on both sides of 


(8.18) we find 
| da 


Ty CE) = a (8.21) 
Eqs. (8.18-21) imply 
2 2 da 
du, a Uy dUy au 1 N, 2 
: ee 
ieee, (IK) gege) + (gyre (L-™) ggg) oO SS (8.22) 

Since 

ou, 

ot be wag See 
we obtain 

2 

du, a uy 1 
auy A : 1/2 a 
(aE, (=e) ge =f (rex ?7? ty 729 (auyae) 2x ax 
an 
i “i783 3/2 2 
a a (1-x) (14x) 797? (au sat)? ax, 

Therefore , (8.21) gives 

1 3 

a - 
Sf (1x) (1-x?) 4(N)2 ax £ 9 


1 - (8.23) 


afi ent barbs 


au 
This proves that the evolution of = is stable in the norm 


du 
(8.14). Finally, the boundedness of =x implies the 
n. 


boundedness of u as will now be show If ou is given 


N’ N 
by (8.19), then 
du N-1 
N a 
as ee 
ox n=0 n n 
where 
(. b -b 
oe eck. ee n+l + 
+ = aca eaaman (n = 1,...,N) 
The boundary condition Uy (-1,t) = 0 requires that 
} 
a = 
v n=1 ” 
Therefore, since i is bounded, so is Un 
x 


In Sec. 9 we present a variety of numerical results for 
the numerical solution of (8.1) by Chebyshev and Legendre spectral 
methods. 


Effect of Boundary Conditions on the Stability of Spectral Methods 


Let us discuss the effect of boundary conditions on the 
stability of the Cheybshev approximations to (8.1). In Sec. 6 it 
was shown that incorrect treatment of the boundary does not affect 
the stability (though it does affect the convergence) of the Fourier- 


Galerkin method. This is not the case for the Chebyshev-spectral 


methods. Let us assume that we solve (8.1) ignoring the boundary 


condition (8.3) and suppose that Uy Oe 0) = Ty Od. The resulting 


21 66- 


system of Galerkin equations for {a} is 


da : 


co - c ) pa (8.24) 


n p=ntl P 


ptn odd 


where a, (0) = ) Eq. (8.24) can easily be solved: a,_, (t) 


nN* 
is a polynomial in t of degree k of the form 


atte ay CS aes. (8.25) 


a This solution is clearly not bounded by any finite power of N. 

: Thus, the Chebyshev methods are algebraically unstable when no 
boundary conditions are applied. 

Ic If we had imposed the boundary condition u(+l,t) = 0 in 
addition to, or instead of, the boundary condition u(-1,t) = 0, 
then Chebyshev-spectral solution to (8.1) would be unstable. 

€ With u(+1,t)=0 instead of (8.3), the Chebyshev-spectral approximations 
to the operator ~9/3x all have eigenvalues with positive real parts 
(that grow as N~*+ ~), Similarly, if we tried to impose the extra 

() boundary condition du(+l,t)/d9x = 0 in addition to u(-1,t)=0 [as 
is frequently done with finite difference methods], an unstable 
scheme would result. 

() The effect of imposing u(+l,t) = 0 in addition to u(-1,t) = 0 
is slightly different for Leqendre-spectral methods. With u(-1,t)= 
u(+1,t)=0, Legendre-spectral methods for solution of (8.1) are 


{ semi-bounded. In fact, 


1 
(v, (LL) v) = =2 J vov/ox ax = 0 
-l 


when v(tl,t) = 0, so these methods are semi-bounded and stable. 


| 
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However, these spectral approximations are not consistent. 


For example, Galerkin approximation involves expansion of 


u(x,t) in terms of the functions bon (x) =Po, (x) = Py (1) 


Songer (X) =Pangy (x) P, (x) that satisfy >, (#1) = 0. 
But 9du/dx cannot, in general, be expanded in terms of the 


functions $,, (x) ; 


The above situations are typical of rapidly converging 
spectral methods. Spectral methods are extremely sensitive to 
the proper formulation of boundary conditions. When proper 
boundary conditions are imposed so the problem is well posed, 
the methods yield very accurate results; when improper boundary 
conditions are mistakenly applied, the methods are likely to be 
explosively unstable. The stability and convergence of spectral 


methods follows very closely that of the exact equations. 


9. Time Differencing 


In previous sections we have investigated the properties of 


spectral approximations to the spatial operator L of the 


differential equation 


In this section we investigate the properties of time-integration 


techniques for the solution of the semi~discrete spectral approx~ 


imations 


Time discretization errors in both finite difference and 


spectral methods are typically much smaller than are spatial 


There are two reasons for this: (i) time 


discretization errors. 


steps are frequently restricted in size by explicit stability 


conditions -- stability of the time integration requires that 


time-differencing errors be small; and (ii) many problems involve 


several space coordinates so any possible efficiency in the 


representation of the spatial variation of the dependent variables 


is quite important to the overall efficiency of the method-- if 


the number of degrees of freedom necessary to describe a certain 


three-dimensional field accurately can be reduced by two in each 


Space direction then the total number of degrees of freedom is 


decreased by a factor 8, but a similar improvement in time 


differencing gives just a factor 2. We will investigate 


here only finite-difference methods of finite~-order accuracy for 
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timewise solution of (9.1) despite the infinite-order 

accuracy in space of many of the spectral methods discussed 

in earlier sections. No efficient, infinite-order accurate 
time-differencing methods for variable coefficient problems - 

are yet Known. The current state-of-the-art of time-integration 
techniques for spectral methods is far from satisfactory on both 
theoretical and practical -grounds and the results to be presented 
here must be regarded as only a .keginning. 

One of our prime goals is to investigate the stability of 
time differencing methods for the solution of (9.1). To do 
this we must first explain how to extend the stability definitions 
given in Sects. 4 and 5. Let uy (x) = @, (x, ndt) be the approx- 
imation to the solution of (10.1) at time nit, where At is a 
time step. Time differencing methods involve approximating 


P ‘ ‘ n+1 
in some way to give a rule for constructing Uy : 


n+l _ n ‘ 
Uy = K, (At) uy ' (9.2) 


where Ky is an operator acting on Uy: Using this rule repetitively 


it follows that 
@ (x,nAt) = (at) 1" 0 
Wy x,ndt) = [Ky J Uy OX ), (9.3) 


where for notational simplicity we assume At fixed. We say that 


(10.2) is strongly stable if 


[| (kK, (At) I™] | < K(nAt) (9.4) | 


u 


for all N and n sufficiently large and At sufficiently small. 

Here K(T) is a finite function of T. We define generalized 

stability by replacing K(T) in (9.3) by n@tPT x (7) as in(5.2). 
A sufficient, though not necessary, condition for strong 


stability (9-4) is ¥ 
[| K, (at) |] -1 < «dt (9.5) 


for some finite « andall At sufficiently small. [If Ky (At) 
is a normal matrix then stability is assured if the eigenvalues 


A of Ky satisfy the von Neumann condition 
max|A| < 1 + Kdt (9.6) 


for sufficiently small At (Richtmyer & Morton 1967). If Ky 
is not normal, then (9-6) is still a necessary, though not 


sufficient, condition for stability in the sense of (9.4). 


The importance of these stability definitions is that they 
lead to the fully discrete form of the equivalence theorem (see 
Sec. 4): a scheme is consistent if 


Ky (At) =< 1 


|(2*——_ - bul + 0 (9.7) 


as N+ and At + 0 for all u ina dense subspace of H; 


scheme is convergent if 


[Jus - u(ndt) || + 0 


as N+ and At +0 for all n_ satisfying 0 < ndt <T and 


; all u(0) eX. The equivalence theorem states that for consistant 


approximations to well-posed problems, stability is equivalent 


to convergence. 


Let us now study the stability properties of some 


specific time-differencing methods. 


Implicit time-integration methods 


Two time-integration methods that are unconditionally 


stable for every algebraically stable spectral method are the 


Crank-Nicolson scheme and the backwards Euler scheme. For any 


semi-discrete spectral approximation (9.1) to uy = Lu, the 


Crank-Nicolson time-differencing scheme is given by 


n+l 


n+l n Uy 


n 
u - u = At Ly ( ) (9.8) 


2 , 


and the backwards Euler scheme is given by 


n+l ni y Atl 
oN ~~ At Ly oN . (9.9) 


To prove that (9.8) or (9.9) is stable, we proceed as follows. 


If (9.1) is algebraically stable there exists a family of 


positive definite Hermitian matrices {H,} such that 


Hy, by + Ly* Hy < o(N) # 


N N N N 


or, equivalently, 


-1 2 
72 yt Hy /? < a(NiT, 


where a(N)<d%nN for some finite d. Substituting 


into (9.8-9), we obtain, respectively, 


v ty nol 


atm, ( S_,4—_), 


Taking the scalar product of (9.10) with 


ny, yntl 


N N°? we get 


v 


My+M,* 


1 n n+1 n 
+¥,, ni (Vi +¥y »)) 


n,,;2. adt n+l, ,2 n 
[17< SSE Live ly 


2 
{| 


(9.12) 


} 


Pe. 


— ere a gg PPS toes 


Therefore, 


1,(;2 
[Iv FA] \f< (1+ 3% gat) ny | 2 
N = poten lity 11") (9.13) 


—_ 


which proves generalized stability for v and, hence, also for 


N 


sg 4, 


WN N N 


Similarly, we may show that the backwards Euler method 


is unconditionally stable. Taking the scalar product of (9.11) 


n+l n 


with Vv + YN we obtain 


N 


so that 


2 


1 n 
[init 


l-aAt 


n+l 
11 Yn 


2 
i 


proving generalized stability of Wyre 
Note that the above proofs show that if a(N) is not a 
1/2 


function of N_ then %; = Hy uy is strongly stable for both 


the Crank- Wicolson and backwards Euler schemes. 


Spectral approximations using Fourier series 


Next, we consider several time integration 


methods for Fourier series spectral approximations to 


with periodic boundary conditions. As shown in Sec. 6, the 


collocation equations are 


du -1 
Cc DCuy (9.15) 


where the matrices 2Nx2N C and D are defined in (6.3). 


The ‘leapfrog' time differencing approximation to (9.15) is 


the explicit two-level scheme 


n+1 n-1l -1 n 
uy - Uy = 2At Cc DCUy (9.16) 
Thee, in the leapfrog scheme 
n n-1 -1 n 
kK (athuy = Uy + 2Atc bcuy ' 


so Ky is a two-level evolution operator since it depends on 
both ee and uN: The definitions of stability, convergence, 


and consistency given above extend easily to this case. 


We shall show that (9.16) is strongly stable provided that 


Ta cammaratateed 


1 


To show this we first recall from Sec. 6 that C is unitary ] 


1 


and D is skew~-Hermitian. Therefore, A = C “pC is also skew- 


Hermitian, and hence normal, so that 


ae ERS cca 


}a}| = 2m(N-1) . 


Now we take the inner product of (9.16) with Uy +u 


to get 


n+l n-l 


+ - 
[u®*2) |? - [fu*][? = 2aerecunt* + ut", aud) , 


; n+ . 5 
since — and uy are real. Since A* = -A, we obtain 


ns gp, Pt ny) )2 _ n+1 
Ul = Muy {{ + hast | 2AtRe (uy 


n 
N ¢ AUy) 


2 “1 -1 
= [lust | + {jak }) = 2AtRe (uy, Au, . 


u )} sg 2° 
N =Uy 


sO uy = Uy - Schwarz' inequality implies that 


n+l Sis n+l 
IRe(un*? ad@y cof fal] Plu®* || | fury 


so that if (9.17) is satisfied, i.e. At||A|| < l-e for some c > 0, 


n+1 n n+l 
|2atRe(uy"”, Aux) | < 2(1-e) Illus dd [Just ; 


Using this result, we obtain 


+1 
ectfay dL + (agli?) + Ged Ciiugt [= [[ub{ 1?) < ut = 02 
(9.18) 


0, 
Since Uy is a bounded function of N (because of the initial 
conditions) , we see that toe is bounded for all N and 
n, proving strong stability. 
Another way to prove that the leapfrog and Crank-Nicolson time 
differencing schemes are strongly stable for (9.15) is to use a } 
modal analysis, which is justified because A is normal. Thus, 


Le uy is an eigenfunction of A with eigenvalue A _, the 


Crank-Nicolson approximation to Ky (At) is 


ae 


0. i sa 0 
K, (At) uy (1 + 5 AAt)/ (1 > AAt) Ww (9.19) 


Since the eigenvalues A of clpe are all pure imaginary, it 
follows that [1K (at) || = 1, so Crank-Nicolson differencing is 
stable. 


Still another time differencing method for solution of (9.15) 
is to use a Runge-Kutta scheme. It easily verified the first and 
second-order Runge-Kutta methods are unstable unless At satisfies 


conditions that are much more restrictive than (9.17). With the 


first-order Euler method 


_ “4 n 


uy + AtAUy , 


stability requires that N7At be bounded as At + 0 [because 


F (At) |J= 1 + 0 (Nn? At?) ] ; with the second-order scheme 
Ky 


antl/2 om pS n 
WwW Uy + 2 AtAu,, 


We ee . - 


n+1 n ~n+1/2 
Uy Uy + AtAu,, r 


stability requires that ni/3 


At be bounded as At + 0. However, 
the third and fourth-order Runge-Kutta methods give conditional 
stability restrictions like (9.17) which we will now derive. 

The third-order Runge-Kutta scheme may be written fof a linear 


equation like (9.}) “as 


Hi ~« ey 2 Sen « n 
a = [I + AtA + 1/2(AtA)” + 1/6(AtA) Ju, = Ky (At) u,, 


(9.20) 


Since K, (At) given by (9,209) is normal, 


| [K (at) | | = max [1 + AAt + 1/2(Adt)* + 1/6(Adt) >| 
A 


where the maximum is taken over all the eigenvalues of A. 


eigenvalues of A are ik with |k| < 2n(N-1), so (9.6) is 


satisfied provided that 


| 
ia 
iq 
| 
7 


tle Ade ao. 


Thus, this method allows time steps that can be Y3 times 

larger than with the leapfrog scheme while maintaining stability. 

However, if the operator A is complicated, the third-order 

Runge-Kutta scheme requires about 3 times as much work as leap- 

from at each time-step, so it is probably not competitive. — 
Similar analysis of the fourth-order Runge-Kutta scheme 

gives the stability condition 


fe (9.22) 


> 


Thus time steps can be nearly three times larger than with 
leapfrog steps. However, fourth-order Runge-Kutta differencing 
requires about four times the work of leapfrog differencing, so 
the scheme is probably not too useful unless very high accuracy 
is desired. 

Now we shall consider time-differencing methods for Fourier 
series spectral approximations to the heat equation with periodic 


boundary conditions: 


aver i aa i 


ace ie a 


The matrix c-* p*c is negative definite. Because (9.19) still 
holds and all eigenvalues A are negative, Crank-Nicolson time 


differencing is unconditionally stable. On the other hand, it 


is easy to show that leapfrog differencing is unconditionally 


unstable. In fact, if uy is an eigenfunction of clp*c with 


eigenvalue .\ < 0 then | 1K, (at) Pal || grows like 


n 
(-AAt + A+ (At) ?) ~e \(ndt) as At + 0 for fixed \ and 


nAt. Since max|A| = 4n7(N-1) 2 grows like N? as N+ ©, 


| 
| 
| 
| 
| 


[1K (at) Puyl | cannot be bounded be a finite function of nAt 
for all N, proving unconditional instability. 

Another way to solve (9.24) is to use a generalized Dufort- 
Frankel scheme 


n+l a n-1 
N 


(9.25) 


If y> a? then this method is unconditionally stable 


(Gottlieb & Gustaffson 1976). 
Similarly, Euler time differencing of (9.24) is conditionally 
stable. Stability requires that At max|\! < 2 or At < (an? (w~1) 277}. 


Higher-order Adams-Bashforth schemes have similar conditional 


stability limits. 


Se 


| 
| 


(| 


Time-differencing for mixed initial-boundary value problems 


Some care is necessary in the formulation of time- 
differencing methods for spectral approximations to mixed 
initial-boundary value problems. The sensitivity of spectral 
methods to the proper formulation of boundary conditions, 
as shown in Sects. 6-8, carries over to the formulation 
of time-differencing methods for these approximations. For 
example, for most mixed initial-boundary value problems leap- 
frog time differencing is unconditionally unstable for spectral 
approximations. Furthermore, explicit time integration methods 
may be unduly restricted by conditional stability requirements 
in spectral approximations. The origin of these severe 
restrictions is the very high resolution of spectral methods 
near boundaries. Thus, it is frequently necessary to combine 
special kinds of implicit time-integration methods with spectral 
approximations in order to maintain high accuracy at reasonable 
computational cost. Several examples will be given later. 

Let us begin by studying time-differencing methods for 
the Chebyshev-spectral approximation to the mixed initial- 


boundary value problem (8.1-3); 


Up 7 “© 0 (~l<x<1l, t>0), (9.27) 
u(x,0) = £(x) (-l<x<l), (9.28) 
u(-1,t) = 0 (t>0). (9.29) 
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In Sec. 8, we proved that various semi-discrete spectral 
approximations to (9.27-29) are algebraically stable. 


Let us first consider the leapfrog time-differencing 


scheme 


n+l _ n-1l n 
Uy =u + 2dt Ly wy, (9.30) 


where wy (X) is the time-discretized approximation to 

Uy (x, ndt), at is the time step, and the semi-discrete 
approximation is Juy/ ot = Ly Uy: 

This scheme is unconditionally unstable for any At as Nr, 
To show this we recall that in Sec. 8. we proved that the 
eigenvalues of Ly have negative real part (see Table 8. 3) 


and that the largest eigenvalue of L has a negative real 


N 
part that grows like N as Nem, Let us rewrite (9.30) > 
in the 2 x 2 block-matrix form 
uns 2 at L I ul 
N N N 
™ (9.31) 
Uy I 0, “* ou 


If the eigenvalues of Ly are denoted aS iy, then the 


eigenvalues of the matrix on the right in 9.31 are 


= uy Abt v 1 + (At) 2 it 


For fixed N and At>0, 


Hy At 

A ae * + Otae*), 

D (9.33) 
~ Thus | 

om ~y ndt 

Es | Men ite (24 afaen) (O<ndt<T, dt+0) 
(9.34) 
Since [Tx cat] papa ty? and there are eigenvalues of 
2 


ten with negative real part of order N“, no inequality of 
the form (9.4) can be satisfied. Thus, leapfrog time 
differencing of the Chebyshev approximations to (9.27-29) 
is unconditionally unstable. 

- There are several conditionally stable explicit time- 
differencing approximations that can be used with spectral 
approximations to (9.27-29). Two examples are the Adams- 


Bashforth scheme 


os uy + ; At Ln - . at -" (9.35) 

and the modified Euler scheme —* 
urtl oul + at Ly (9. 36a) 
=" = un + $ at Ly, un + 5 at _e (9.36b) 


The modified Euler scheme (9.36) is in practice stable provided 


the stability condition 


8 
At < 
Nn? 


is satisfied. A similar stability condition holds for the 


Adams-Bashforth scheme. 


(9.37) depends 


The fact that the stability limit in 
2 


on 1/N rather than 1/N is not very surprising 


because the Chebyshev collocation points 
are spaced by a distance of order isn? near the boundaries. 


cos tn/N 


Since the wave speed in (9.27) is 1 the wave propagates from 


one grid point to the next in a time of order 1/N? so time 


steps must be smaller than this to maintain explicit stability. 


The explicit stability restriction (9.37) for Chebyshev- 


spectral methods with N polynomials should be contrasted with 


the corresponding stability conditions for finite difference 


approximations to (9.27-29), with N gridpoints uniformly spaced 


in the interval -l<x<l, the grid spacing is 2/N_ so the 
Courant stability condition is At <2/N. As N+, this 


stability condition on finite difference schemes is much 


weaker than the condition (9.37) on the spectral approximations. 


A semi-implicit technique that permits stable time-differencing 


with spectral methods with a stability condition like that 


of finite-difference schemes will be discussed later in this 


section. 


In order to prove that the modified Euler method (9.36) 
is stable, we begin by noting that (9.36) is equivalent to the 


second-order Taylor series method 


uptl = (7 + at Ly + Z(At)? 147) uy = Gy uN (9.38) 


A sufficient condition for algebraic stability of (9.38) is 


the existence of positive-definite symmetric matrices S, 


such that 
T (9.39a) 
Gy Sy Sy < Sy 
and the condition number of Sy Satisfies 
Hs IT Istt)) = o¢n8) (yee) 
N N ; (9.39b) 
for some finite 8 . If (9.39) holds then 
T,n n T,n-1 n-l 
or 
-1/2 42 /2 n .-1/2 ———— 
Sy agra ere gine 
Therefore, ; 


1/2 nN o-1/2 
ey we?) | <1. 


PRISE We aie 


so that 
0 -1/2 1/2 n .~1/2 
[uBiy = [feoyy® ugtistisy/21 tL isy/ 7a)” sy? 


Ls 771TH ag tt = 0 xP} fugit) are). 


To complete the stability proof we must investigate 
under what conditions matrices Sy Satisfying (9.39) exist. 


One choice for S, is just the Liapounov matrices of Lyi 


these matrices satisfy 
jin+ i c= +2 (9.40) 
na * “uy Sn 


It was shown in Sec.8 that the Liapounov matrices for spectral 


approximations to (9.27-29) have algebraically bounded condition 


number. Using (9.38) , we obtain 


1 


T T 21,2 1 2 2 
Gi Sy Gy = [T+ At Ly + 5 (Ot) “(Ly ) Sq [T+AtLyty (At) ~ (Ly) ] 


or 


T j T 
Gy Sy Sy = Sy + At (Ly Sy + Sy Ly) 


—— oe 20 T 2 
- x (At) ( (Ly) Sy + 2h, Sx Ly + Sy Ly ] 


TS 
+2 (aty?tingy” On ON + LY 


2 1 4,.a,¢ 2 
Sy byl + g (at) (Ly) Sy by - 


From (9.40), it follows that 


aoe 


a 


2 
Sy AtI - (At) (Ly 


3 ,.T 1 4 Pes Zz 
ba Ly at z (At) (Ly) S.. L 


1 
= & tee) NN 


Thus, (9.39a) is satisfied provided that 


T 
-At (Ly + Ly) <2 (9.41) 


(9.4 2) 


If (9.41-42) are satisfied then the modified Fuler method for 
(9.27-29) is algebraically stable. 

At first , it may appear that the stability condition 
(9.42) is much more severe than the stability condition 


(9.44). In fact, we know from Sec. 8 that 


- 2 ‘ < 
In, 11 = own), IIs, 1] = 0@) (N+) , 


so that (9.42) seems to require that At=0(1/N*) as 
Neo, However, the stability condition (9.42) is no 
more restrictive than the stability condition (9.41) 
To see this we use (9.40) written in the form 

1 Le eB 


+ (Ly) L SS. b= = ft 


. = 
De NS Ra an N ON UN 


N “NN N 


to obtain the representation [see (5.13)] 


expl(Lo")"t] exp[L, t] ae. 
(9.43) 


It may be shown that the norm of the integrand of (9.43) 
is O(1) as N+ ~ for t = 0 (N2) and that the norm decays 


rapidly to zero as t + », Therefore, 


T 
ty Sy Eyl] = 9 OV) are) (9.44) 


showing that the stability condition (9.42) is of the form 


At = O(1/N?). 
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semi-implicit methods 


When explicit time-stepping methods are used to solve 


semi-discrete spectral equations for the hyperbolic problem 


du gu _ Sh 
xed a(x) a 0 (-l<x<1) (9.45) 
with appropriate boundary conditions [that depend on the sign 


of a(x)], there result stability conditions of the form 


1 1 


See ee ae ee 
N*|a(1) | N“|a(-1) | 


i 
Nmaxyatsy|? (9-46) 


| x| <1 


These stability limits can be derived heuristically from the 


Courant stability condition 


Ax 
= (9.47) 


t 
O° & Teese! 


e 


where ere is the effective wave propagation speed in a 
direction in which there is effective grid resolution AX ¢¢- 
Near the boundaries x=+l, the Chebyshev-spectral methods have 


a 2 co i = e 7 
resolution AX o¢¢ O(1/N*) as N-> while aore a(tl); in 


the interior of -1l<x<l, Chebyshev series have effective resolution 


AX.e¢ = 04 x ) as N+» while the largest wave spped is max|a(x)|. 


Thus, (9.47) implies (9.46) for the Chebyshev-spectral methods. 


The stability condition (9.46) is too severe for many 


applications because it requires that At = 0(1/N?). 


In order to relax this severe constraint, we use a semi-implicit 
method in which the propagation through the high-resolution boundary 
is treated implicitly, but the propagation through the interior is 
treated explicitly. 

One possible semi-implicit scheme is the following two-step 
Method. Let Ly be the Chebyshev-spectral approximation to 
“a (x) sk with appropriate boundary conditions applied, and 
Li Ly be the Chebyshev spectral approximations to the 
constant coefficient wave operators - a(+l)9/92x, -a(-1)9/9x, 


respectively, again with appropriate boundary conditions applied. 


A semi-implicit two-step scheme is given by 


A 1 
nty - it m3 -, on 
Uy? Athy uy = us gat (Ly ~ Ly) uy (9.48a) 
1 1 
uptt - dat ue ubtt = utter + dat (L - f) unte —— (9.48b) 


The scheme (9.48) is stable if the stability condition 


at < WRax Tate] nes 


is satisfied. 
The condition (9.49) is sufficient to ensure stability, 


but the semi-implicit scheme (9.48) may be stable even if 


(9.49) is violated. If maxja(x)|<|a(l)] or max{a(x)|<{a(-1)|, 
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(9.48) is usually unconditionally stable for sufficiently 
large N (see Sec. 8 of Orszag 1974). The implementation 
of (9.48) on a computer is straightforward and efficient; the 
properties of Chebyshev polynomials summarized in the Appendix 
show that the implicit equations (9.48) are essentially tridiagonal 
matrix equations. 

The reason that the semi-implicit method outlined above 
does not have a stability restriction like At = 0 (1/N?) can be 
understood as follows. By subtracting Ly and Ly in succeeding 


half-time-steps, the explicit part of the calculation is similar 


to that in solving an equation of the form 


as (1-x7) b (x) au nD (9.50) 


where the wave speed vanishes at x=+l. If b(x) =b, a 


constant, the Chebyshev-tau equations for (9.50) are just 


a 1 = 
— = 2 c, b [(n-1) ain-1| (n+l) a (9.51) 


n+l ] 


where Cy = 2 and * 1 for n>0. By Gerschgorin's theorem, 


| | for (9.51) satisfies 


fi 


[|L.|| = O(bN) = (Neo), (9.52) 


I 


so the explicit time step restriction is At = 0(1/bN) as 


N?o 


We note that Chebyshev-spectral approximations to (9.50) 
are stable when no boundary conditions are applied. In fact, 
using Galerkin approximation and the Chebyshev inner product, 


we obtain 


Oy 
— } om -@-s 


au, 
N Pt 
(Uy. << + b(1l-x”) 


Therefore, 


Iu ted 12 < el PIE)juyco 11? 


Proving stability. 

There are other attractive semi-implicit schemes for (9.45). 
For example, suppose a(x) is one-signed, say a(x)>0, and let 
aiax * maxa(x) . Define i. as the Chebyshev approximation to 
- 4 ee wx with boundary conditions imposed at x =-l. A 


semi-implicit Chebyshev spectral scheme for (9.45) is 


n+l max n+l _ ,n = 
Uy «At Ly Uy uy + At (Ly 


The scheme (9.53) is usually unconditionally stable and avoids 


the severe time step restriction (9.46). It is also easy to 


max 


implement efficiently because Ly is a Chebyshev approximation to 


a constant-coefficient wave operator? 
The same kind of trick stabilizes spectral methods for non- 


linear equations. For example, if we are solving the equation 


during a time interval in which u(x,t) is smooth (no shock waves), 


then we may use the semi-implicit scheme 


in which the terms on the left are treated implicitly in time, 


while those on the right are treated explicitly. Here ee is 


an estimate of the largest value of u(x,t). Similar semi-implicit 
methods are effective in eliminating (or at least relaxing) time- 
step restrictions for finite-difference methods. The key idea 

is to recognize the source term of a numerical instability and then 
to approximate it by a simple expression that can easily be treated 
implicitly. 

ied ae 


t , ; 
‘e Haidvogel has pointed out that the semi-implicit scheme 
-53) with L replaced by a Chebyshev spectral approximation 
to &(bxt+c)3/3X, where bte = a(+l), c-b = a(-1), is also 
Stable under the weak restriction (9.49). The resulting 


implicit equations are still tridiagonal [see (A.9), (A.18)]. 
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Several other examples of semi-implicit methods should make 
the general technique clear. For the variable coefficient heat 


equation 


= k (x) ws (-1<¢xX<1) 
with suitable boundary conditions at x =!l and k(x)>0, 
Chebyshev-spectral methods give explicit time-step stability 
conditions of the form 
P 1 
At < min (otic ne erat. Soka =?) 


’ 
kK(-2)N4 ok) N4 Os N* max k(x) 


|x] <2 


: pra 4 
The very severe time step restriction that At = O(1/N ) as 
Nee is due to the high resolution of Chebyshev series near the 
boundaries x = +1. To avoid this problem we can use a semi-implicit 


method. Let Ly be the Chebyshev-spectral approximation to 


k(x) 92/ax" and let a be the Chebyshev-spectral approximation 


to dy 97 /9%7 where Pine = maxk(x) . The semi-~implicit 


2 max 
max 


scheme (9.53) with Ly defined in this way 


seems to be unconditionally stable (Orszag 1974) and certainly 

does not have any stability restrictions of the form (9.54). 
Finally, we comment on the need for implicit or semi-implicit 

schemes in multi-dimensional problems. If we wish to solve the 


Navier-Stokes equations 


Va hae Sela w e 


for incompressible fluid flow, the treatment of the various 
terms should guided closely by the type of stability restrictions 
they impose. 

If v= 0 then we need only consider the types of stability 
restrictions induced by the advective term -U.Vu and by the 
pressure term -Vp; we will not discuss the effect of the 
pressure because it is closely connected to the incompressiblity 
condition V.u=0 and is not relevant to the semi-implicit ideas 
discussed here. At a boundary of the flow, it is appropriate to 
specify boundary conditions on u.n where fh is the normal 
to the boundary. If the boundary is solid and stationary, then 
u.n=0 and we are in a situation similar to that modelled by 
(9.50). The effective convective speed normal to the boundary 
vanishes, so spectral methods exhibit no unusual time stepping 
restrictions. However, if fluid is being blown into or sucked 
out of the boundary so u.n# 0, then semi~implicit methods must 
be applied to avoid unreasonably restrictive conditions like 
(9.46) on the time steps. 

If v>0, then the viscous terms in the Navier-Stokes equations 
should be treated implicitly to avoid unreasonable time step 
restrictions due to the high resolution of spectral approximations 


near the boundary. 


-195- 


i 
‘| 
] 


10. Efficient Implementation of Spectral Methods 


There are two aspects of the efficient implementation of 
spectral methods that we discuss here: (i) evaluation of 
derivatives; (ii) evaluation of nonlinear and nonconstant 
coefficient terms; (iii) roundoff errors. More details on 


these matters are given elsewhere (see the References) . 


Evaluation of derivatives 

An efficient procedure to obtain the expansion coefficients 
of derivatives of a function f(x) in terms of the expansion 
coefficients of f(x) is to use recurrence relations. For 


example, to evaluate the term 


} 

Ss, = pa 

n p=ntl Pp 
pt+n odd 


that appears in the Chebyshev equations (2.11), (2.19), and 


(2.32),we use the recurrence 


ji nt2 * (Mtlanyy (O<n<N-1) (10.1) 


with Sy = Snel = 0. In this way, §. is evaluated for all n 


using only N arithmetic operations. The existence of the recurrence 
relation (10.1) is ensured by the recurrence property 
' 
Th+1 n- 


Sa 


r a ar 


fee 


t 


satisfied by the Chebyshev polynomials. Similarly, it is possible 


to derive recurrence relations to evaluate efficiently the 
coefficients of arbitrary derivatives of functions expanded in 
Chebyshev and other classical polynomial expansions. 

Evaluation of nonlinear and nonconstant coefficeint terms 

The most efficient way to evaluate nonlinear and general 
nonconstant terms in spectral approximations is to apply transform 
methods. The key idea is to apply fast Fourier ea Med and other 
transforms to transform efficiently between spectral representations 
of a function f(x) and physical-space representations of f(x). 
With Chebyshev polynomial expansions, fast Fourier transforms permit 
the evaluation of arbitrary nonlinear and nonconstant coefficients 
terms in order N log N arithmetic operations. 

In general, collocation methods give approximations to nonlinear 
and nonconstant coefficient problems that can be more efficiently 
implemented than Galerkin or tau approximations. Collocation is 
recommended for these problems. For example, the solution of the 


hyperbolic problem 


+ 
Sot eee ee em Bie th (-L exe 4, ts OP, 020.3) 
t ox - ~ 
u(-1,t) = 0, 


would be difficult using Galerkin or tau approximation but is 


straightforward using collocation methods. 


Let us explain how to march the solution to (10.2) forward 
by one time step efficiently using Chebyshev collocation. We introduc: 


the Ntl collocation points x5 = cos 1j/N (j = 0,...,N) 


and represent the current solution 5 as 
N an 
u = z a, cos L2. 2 (10723) | 


Then we invert (10.3) by the fast Fourier transform to 


obtain ay for n= 0,1...-.,N and calculate 
a = 28/¢, 
by (10.1). Next we evaluate 


N ‘ 
au « ££ ae cos 74 (10.4) 


using the fast Fourier transform. Finally, we evaluate 
exp(u5+x5) (9u/2x) 5 at each of the 'grid' points x5 


and use the results to march the solution forward to the next 


time step. 


For quadratically nonlinear differential equations, like the 
Navier-Stokes equations of incompressible fluid dynamics, Galerkin 


and tau approximations are workable but normally require at least 


twice the computational work of collocation approximation. However, 
Galerkin approximation is sometimes very attractive because it gives 
approximations that are conservative and have no so-called aliasing 
errors (see Orszag 197lc, 1972 for a more complete discussion of 
these properties). 

Roundoff Errors 

Transform methods normally give no appreciable amplification 
of roundoff errors. In fact, the evaluation of convolution-like 
sums using fast Fourier transforms often gives results with much 
smaller roundoff error than would be obtained if the convolution 
sums were evaluated directly. 

On the other hand, the use of recurrence relations to evaluate 
derivatives can sometimes be a source of large roundoff errors. 
In this case , it is often best to convert the problem being solved 
into a new one that is numerically well-conditioned. An example of 


such a transformation is given below. 


Example 10.1: Solution of y"-ky=f(x) by Chebyshev polynomials 
The boundary-value problem 


y" SRY © Ste ~1<x<l (10.2) 
y(-1) = A, y(l) =B 


can be solved using a Chebyshev-tau approximation. The 
resulting approximation Y yy 6) is given by (see Appendix) 


N 
Yq (*) ts ay 7, (x) (10.3) 
f 
| 
| 
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x 2/42 
t posh i) 4. ca fh (O<n<N-2) (10.4) 


e, p=n+2 Pp 
ptn even 
N rm N 
t @iy" &. = B, ») a =B, (10.5) 
n=0 n n=0 n 


where {£1} are the Chebyshev series coefficients of £ (x). 
The solution of the system (10.4-5) for the Chebyshev 


coefficients {a,} may be done in several ways. One obvious 
way to do this efficiently is to write 


ee. + aap + Ba, ANE (10.6) 
(1) err i | all) ws: 
Here ap satisfies a, = ay) = 0 and 
1 OY sald) 
+ }  plp*-n ya -ka) = £ (0c n<N-2), 
n p=n+2 n n 
’ (2) nen (2) (2) 
while a, satisfies a, = 1, a, = 0 and 
1 N 2 (2) (2) 
. ee | p(p°-n ya, ~ ka, = 0 (O<n<N-2), 
fy p=n+2 , 
3 es 3 3 
and a! ) satisfies an ) = 0, al?) = 1, and 
1. 4 i2) a (8) 
a ) p(p- -n wed ka = 0 (0<n<N-2). 
nN p=n+2 
ptn even 


Each of the solutions "al ae ‘ ai » May be found 


using roughly N operations by backwards recurrence. When 
the constants a and 8 in (10.6) are chosen so that the 
boundary conditions (10.5) are satisfied, an given by 
(10.6) satisfies (10.4-5). 

The above procedure is efficient but it is not usually 
numerically stable. Roundoff errors multiply rapidly so that 
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a, may have little significance. 

A better procedure is to first convert (10.4-5) into a 
nearly tridiagonal system of equations. It may be shown that 
(10.4-5) is equivalent to the system 


—na=2 n+2 n+4 
S46 - () ern a_+ a 
4n(n-1) ® 2 2(n“-1) n ~ 4n(ntI) “n+2 
c £ e £ e £ 
ao _o72 n-2 . —nt2 on n+4_n+2 
~ 4n(n-1) Bie AE + -aery  (2sns) (10.7) 


with the boundary conditions (10.5) still applied. Here 

Cy=2, c,=1 for n>0O and e,=1 for n<N, e,=0 for n>N. The 

system (10.5), (10.7) may be solved by standard banded matrix 
techniques in roughly the number of operations required to 

solve pentadiagonal systems of equations. The equations in the 
form (10.7) are essentially diagonally dominant so no appreciable 
accumulation of roundoff errors occurs. This technique for 
solution of (10.2) is very useful in implementing implicit 
spectral methods for dissipative terms and for solving Poisson- 
like equations (see Sec. 14). 
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| ll. Numerical Results for Hyperbolic Problems 


We begin by presenting numerical results for spectral approximations 


to the problem 


u +tu. = 0 (-1<x<1,t>0) (11.1) 
t x meses 


u(x,0) = 0 ,u(-1,t) = g(t), (11.2) 
whose exact solution is 


u(x,t) = 
(x>t-1), (11.3) 


a =x - 1) (x<t-1) 
If g(t) is smooth, u(x,t) is smooth for [x|<1 when t>2; when 
t<2, u(x,t) is not smooth at x=t-l. 

In Sec. 2. we explained how to obtain semi-discrete Galerkin, 
tau, and collocation approximation to (11.1-2) using either 
Chebyshev or Legendre polynomial expansions. In Sec. 9, we showed 
that either Adams-Bashforth or modified Euler time differencing gives 
stable and convergent results for these spectral approximations. The 
numerical results cited in this Section were obtained by Adams-Bashforth 
time-differencing; time steps were chosen small enough that time- 
differencing errors are negligible. 


Comparison of Chebyshev and Legendre Polynomial Spectral Methods for 


Smooth Solutions 


When g(t) =-sinMnmt, the solution (11.3) has M_ complete 
waves within |x|<1 when t>2. As argued in Sec. 3, we expect that 
accurate results will be obtained only if N>M" polynomials are 


retained. 


Ga 


In Fig. 11.1, we plot the root-mean-square error for |xj|<1 
averaged in time for 4<t<4.4 obtained using the Chebyshev 
approximations to (11.1-2) when g(t) = -sin5mt. In this time 
interval, u(x,t) is smooth for |x] <1. Observe that the errors 
decrease exponentially fast when N>5r. Also observe that when the 
spectral approximations are relatively inaccurate (errors greater 
than roughly 10%), Galerkin approximation is most accurate followed 
by collocation and then tau. On the other hand, when the spectral 
approximations are very accurate (errors less than roughly 0.58), 
tau approximation is most accurate followed by Galerkin and 
collocation. This behavior seems typical. Also observe from 
Fig. 11.1 that all three spectral approximations are nearly as 
accurate as the best (rms) Chebyshev approximation; in fact, tau 
approximation with N+l polynomials is usually more accurate than 
the best approximation with N polynomials. Here the best (rms) 
Chebyshev Seoeietye is that Nth degree polynomial that 


ei 


minimizes f |u - ul? (1-x) 
-1 


In Fig. 11.2, we make similar comparisons of the error in 
spectral approximations using Legendre series for the problem 
(10.1-2) with g(t) = -sin5™t. Here too the errors decrease 
exponentially fast when N257. Again, tau approximation is more 
accurate than Galerkin when both are very accurate, while it is 
less accurate when both are relatively inaccurate. Also, tau 
approximation with N+l polynomials and Galerkin approximations 
with N+2 polynomials are more accurate than the best Legendre 
approximation with N polynomials. Here the best Legendre 


approximation is that Nth degree polynomial that minimizes 


2 
iB Ju, - ul “dx. 
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that the errors decrease rapidly for N > 5m, 
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Fig. 11.1. A plot of the Lo-errors in Chebyshev-spectral 
solution of (11.1-2) with g(t) = -sin 51t. The errors 

are averaged in time over the interval 4<t<4.4; the exact 
solution u(x,t) = sin 5(x+l-t) is smooth throughout this 
time interval. The best (rms) approximation is given by 
(3.41) with M = 5, a = l=-t truncated after Ty (Xx) - Observe 
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Fig. 11.2. Same as Fig. 11.1 except for Legendre-spectral 
solution of (11.1-2) with g(t) = -sin 5mt. Here the best 

(rms) approximation is given by Aa 45) with M = 5, a = l-t 
truncated after Py (x). 
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In Fig. 11.3-4 we plot the error Ey (Xt) in the best 
Chebyshev polynomial approximation to sin5%(x+l-t) at t=4. 
Observe that Ey (Xt) is nearly an ‘equal ripple' approximation 
(Acton 1970) so uy (x,t) is nearly a minimax approximation. 

In Figs. 11.5-8 we plot the errors Ey (Xt) versus x 
at t=4 obtained by numerical solution of Chebyshev spectral 
approximations to (11.1-2). As N increases, the tau method 
gives the cissest approximation to an equal-ripple error, which 
is consistent with the result shown in Fig. 11.1 that tau approximation 
gives the smallest errors at high accuracy. 

In Figs. 11.9-10, we plot the error in the best 
Legendre polynomial approximation to sin51(x +1l-t) at t=4. 
Observe that E yy (Ret) has large errors near the boundaries 
x =+t1l. By comparing the results plotted in Figs. 11.3-4 with 
those plotted in Figs. 11.9-10, we conclude that the best Chebyshev 
polynomial approximation is closer to an equal ripple approximation 
than is the best Legendre polynomial approximation. Even though the 
best Legendre polynomial approximation to u(x,t) gives the smallest 
mean-square error to u, the best Chebyshev polynomial approximation 
usually gives a smaller value of the maximum pointwise (L,) error. 
The large errors of the best Legendre approximation are concentrated 
near the boundaries x=+l, while the Chebyshev weight function 


-1/2 


(1-27) tends to distribute the errors in the best Chebyshev 


approximation uniformly throughout -1l<x<l. 
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In Figs. 11.11-13, we plot the errors € (xt) at t=4 
obtained by numerical solution of Legendre spectral approximations 
to (11.1-2). As for Chebyshev-spectral approximations, the error 
in Legendre-tau approximation is smaller than that in Legendre- 
Galerkin approximation. 

One important feature of Legendre-spectral approximation is 
that the spatial distribution of the error in tau and Galerkin 
approximation plotted in Figs. 11.11-13 differs markedly from 
the spatial distribution of the error in the best Legendre 
polynomial approximations plotted in Figs. 11.9-10. The boundary 
errors in the best Ly approximation are relatively large while 
the boundary errors are relatively smaller in the spectral 
approximations. 

The boundary (endpoint) errors in Legendre-tau approximation 
exhibit ‘superconvergence’ in the sense that they go to zero much 
faster than either the Ly - errors or the Lo and endpoint errors 
of Chebyshev-tau approximation. This fact is illustrated in Fig.11.14 
where we plot the Lo and endpoint errors of Legendre-tau and 
Chebyshev-tau spectral approximations to the solution of (11.1-2) 
with g(t) =-sin5nt. Here the endpoint error is | u, (+1,t)-u(+1,t)| 
at the outflow boundary x= +1. 

Several features of the results plotted in Fig. 11.14 deserve 
comment. First, although the maximum error of the best N-term 
Chebyshev polynomial approximation is smaller than the maximum 
error of the best Legendre polynomial approximation to u(x,t) by 
roughly a factor 1//YN [see (3.38) and (3.39)], the maximum error 


of the Legendre-tau approximation is smaller than the maximum error 
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Fig. 11.14 A comparison of the Chebyshev-tau and Legendre- 
and endpoint (x = +1) errors for the solution to (11.1-2) 
(t} = -sin 5rt. 
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of the Chebyshev-tau approximation. Second, the endpoint error at 
x=l1 of the Legendre-tau approximation goes to zero like the square 
of the endpoint error of the Chebyshev-tau approximation. This 
remarkable behavior of endpoint errors in Legendre-polynomial 
approximations was found originally by Lanczos in a slightly 
different context [Lanczos 1966 (p. 156), 1973]. 

A mathematical analysis of the errors of spectral approximations 
to (11.1-2) has been given recently by Dubiner (1977). Dubiner's 
results include: (a) asymptotic estimates of the errors incurred 
by the various spectral methods, including error oscillations when 
the solution is smooth; (b) a complete boundary layer description 
of the decay of large errors due to discontinuities after the 
discontinuities propagate out of the computational domain; (c) 
analysis of the behavior of the tau-function t(t) in (2.34). Dubiner 
has analyzed a variety of spectral methods for (11.1-2) based on 
expansions in general Jacobi polynomials. His ingenious analyses 
of tau methods should permit more complete analysis of these 
methods than possible using earlier a posteriori analysis (see Fox 
& Parker 1968 for examples of a posteriori error analysis of tau 
methods) . 

Mesh Refinement 

Sometimes it is useful to split up a domain into several 
subdomains and then use numerical methods of different spatial 
resolution in each. For example, in limited-area numerical weather 
forecasting near a metropolitan area, it may be desirable to have 
much finer resolution in a small region than is practical globally. 


One way to do this is to solve the problem separately on each 
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of several subdomains and then to match the numerical solutions 
so obtained across subdomain boundaries. As a model of this 


Procedure we consider the problem 


+u. = 0 (-l<x<l, t>0) (11.4a) 


u(-1,t) = g(t), (11.4b) 


v. + v= 0 (1<x<3,t>0) (11.5a) 
v(l+,t) = u(l-,t). (11. 5b) 
with finite difference methods, the accurate solution of the coupled 


system (11.4.5) using different grids for -l<x<l_ than for 


1<x<3 may be troublesome. Inaccurate results or even numerical 


instabilities can result from the matching (Browning, Kreiss & Oliger 


1973). Because grids with different grid separations have different 
dispersion characteristics for waves propagating on the grid, 
waves can reflect from the boundary at x=l1 and cause large 
errors. 

Spectral methods are attractive for the solution of mesh 
refinement problems like (11.4-5) because they give small endpoint 
errors. For example, the Chebyshev-tau approximation to (11.4-5) 
with N+l polynomials to represent the solution for -1l<x<l and 
M+l polynomials for 1<x<3 is given by 


N 


(x,t) = JF a(t) T(x) (-1<x<1) (11.6) 
oN 9 2 n 
n= 
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M 
vy(x.t) = J balt) F,(x-2) (1<x< 3) (11.7) 
m=0 
N 
act EOE 5 y pa (0<n<N-1) (11.8) 
dt on p=n+l P 
ptn odd 
N 
a gare saa (O<meM-1) (11.9) 
dt Cn p=m+1l P 
p+m odd 
N n 
y (-1)2 a = g(t) (11.10) 
n=0 
M m N 
er ee (11.11) 
m= n= 


It may easily be shown that if g(t) is smooth, the solution to 
(11.6-11) converges to the solution of (11.4-5) throughout 
-l<x<3 faster than any finite power 1/N or 1/M as N, M+, 


The solutions for -l<x<l and 1<x<3 match without the 
necessity of imposing any matching conditions in addition to 


(11.5b). Because no spurious downstream boundary conditions 


are applied at x=+l on the wave propagating in the interval 


-1<x<l, there are no reflected waves. 


One more example of a refined mesh spectral calculation 


, is instructive. Consider the heat equation problem 
tu. au 
7 -1<x<1 (11.12a) 
; ax 
j | 
] 
ov a*y 2 
= 1<x<3 (11.12b) 
f 
u(-1,t) = £(t), v(3,t)=g(t) (11.12c) 
' u(lyt) = v( 1+,t), ou Q-,t)= 2 (+,t) (11.124) 


where (11.12a) follows by requiring continuity of temperature and 
. heat flux across the boundary at x=l. A Chebyshev~-tau approximation 


to (11.12) is given by 


N 
P u(x,t) = JF ap (t)T (x) (-1<x<1) (11.13a) 
n=0 us 
M 
v (x,t) * be b(t) T (x-2) (1<x< 3) (11.13b) 
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da N 


poe f p(p’-n*) a (Ocn<N-2) (11.13c) 
_— eae p 
p+n even 
ab M oe 
re = = aes p(p -m )b,, (O<p<M-2) (11.134) 
m ptm odd 
N n M 
y (lp a, = lt), J b= gilt) (11.13e) 
n=0 m=0 
N M N M 
Yo oa= 2 (20. 2 n’a = - 2 (-2)"n%b,. 
n=0 7 m=0 a pg m=0 m 


(11.13£) 


It may be shown as in Example 7.1(v) that this approximation is 


semi-bounded and hence stable and convergent. 


Discontinuities 

When t<2, the solution (11.3) to (11.1-2) is not 
smooth at x=t-l; if g(t) = sinMnt, the solution has a 
discontinuous derivative. This discontinuity seriously degrades 
the rate of convergence of spectral approximations near the 
discontinuity. Nevertheless, spectral approximations are still 
normally much more accurate than finite-difference approximations 
to the same problem. Orszag & Jayne (1974) give comparisons 
between finite-difference and spectral approximations to 
discontinuous solutions; in particular, they argue that if the 
pth derivative of the solution is discontinuous, the rate of 
convergence of Chebyshev-spectral approximations to (11.1-3) for 
t<2 is of order 1/NP as N-**. Dubiner (1977) has given a detailed 


asymptotic analysis of this problem. His results include detailed 


behavior of the error for all x and t and are in good agreement 
with numerical solutions. 

One of the attractive features of spectral methods for problems 
with discontinuities is that the region of large errors is more 
localized near the discontinuity than in finite-difference methods. 
Thus, it should be possible to eliminate oscillations near the 
discontinuity using less dissipation than is required when finite 
difference methods are used. A comparison of the error in Chebyshev- 
tau and second and fourth-order solutions of (11.1-2) for t<2 is 
given in Fig. 11.15. 

Another interesting way to use spectral methods for problems 
with discontinuous solutions has been suggested by Boris & Book 
(1976). The “optimal flux-corrected transport" approximation 
gives good resolution of discontinuities without introduction of 
unphysical numerical oscillations near the discontinuity. The idea 
is to add in an artificial diffusion to smooth the discontinuity and 
then to ‘anti-diffuse' the resulting solution in such a way that no 


new oscillations or maxima/minima are produced. 


Comparison with Finite Difference Methods 


Finite-difference approximations to (11.1-2) must be 
formulated carefully near the boundaries x =+ 1. For example, 


| the fourth-order semi-discrete approximation 


- + 


du 8(u -u )-u u 
, “ #1 j-1 on te 4 
12Ax 


where u,(t) = u(jAx,t), must be modified at x=-1l+Ax, 1-Ax,1l 


j 
because u(-l-Ax,t), u, +Ax,t), u(1+2Ax,t) all lie outside the 


Emax /m 


~h 2/3 
second -order 


~h 4/5 
fourth -order 


~h 
spectral 


2 S59 10 20 50 1l00 200 500 I000 


Fig, 11.15 Maximum pointwise errors in the solution of (11.1-2) 
with g{t) = sin mrt at t = 1.2 when the discontinuity in the 
Solution (11.3) is at x = 0.2. Here e€ is the maximum error. 


max 
and N = 2/h is the number of grid points, or of Chebyshev polynomials 


in the spectral method. 
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Computational domain -1l<x<l. Kreiss & Oliger (1973) discuss 


methods to formulate difference approximations at these grid 
points. However, it is not known how to formulate appropriate 
‘boundary' conditions for arbitrary order difference schemes. 
This difficulty is an artifact of difference methods; a fourth- 
order difference quation requires 4'boundary' conditions while 
only 1 condition (11.2) is properly imposed on the first-order 
differential equation (11.1). 

On the other hand, properly eorauketed spAceral methods 
require no 'spurrous' boundary conditions. Indeed, the imposition 
of a spurious boundary condition on a spectral approximation to 
(11.1), like du/ax = 0 at xX =+1, will induce an unconditional 
instability (see Sects. 8,12). The mathematics of spectral 
approximations follows closely the mathematics of the differential 
equation being solved. 

Spectral approximations also require considerably fewer degrees 
of freedom to achieve accurate results than are required by 
difference methods. A comparison for the problem (11.1-2) is 
given in Table 11.1 for late times at which the solution is smooth. 

In Figs. 11.16-19 we show three-dimensional perspective plots 
of the solution to a simple two-dimensional hyperbolic problem with 


periodic boundary conditions 


dA(X,y,t) _ dA (x,y,t) dA(x,y,t) si 
BAUErYet) _ y BACE YE) 4 x SACYE) 2 9 (11.14) 
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Table 11.1 


Chebyshev-tau 
M 


Second-order 
N M € 


€ 
rs) 


(rms) errors for the solution of (11.1-2) with 


Table 11.1. Ly 
g(t) = sinMmt. The errors listed are measured at t=5 when the 
solution (11.3) is smooth. Time differencing errors are negligible 
and N is the number of grid points or Chebyshev polynomials. 
Observe that to achieve a 1% error, the second-order method requires 


N/M240, the fourth-order method requires N/M215, while the 


Chebyshev-tau method requires N/Mpr. 


with 


A(xt2u,y+27,t) = A(x,y,t). 


The solution to (11.14) is constant along the characteristics 
x+iy= (Xo + iy,)ett. Therefore, A(x,y,2n7)= A(x,y,0) so 
the solution should keep A unchanged after a time 27. In 
Fig. 11.16, we plot the initial conditions used for the calculation 
whose results are plotted in Figs. 11.17-19. In Fig. 11.17 we plot 
the results at t=2n of a second-order centered space difference 
scheme; in Fig. 11.18 we plot the results of a fourth-order scheme 
and in Fig. 11.19 we plot the results of a spectral calculation 
using the Fourier expansion 

A(x,y,t) = J y  atk,p,tyet***4Py | 

|k P 

All three calculations used the same number of degrees of freedom 
but the differences in accuracy are striking. The Fourier-spectral 
method works well even though the convecting velocity (-y,x) in 
(11.14) has jump discontinuities at x=+2n, y=+21. The dominant 
error in all three calculations originates from the 'corners' of 


the initial A(x,y,0) distribution; thus error appears as a large 


lagging phase error in the finite difference solutions which explanins 


the 'wakes' of large errors following the remnants of A(x,y,2T). 


Higher-Order Wave Equations 


The mixed initial-boundary value 


2 2 
292 a  (-1<x<1,t>0) (11.15) 
at ox 


KG 
<> 


Fig. 11.16 A perspective plot of the A(x,y,t=0) used in a 
numerical test of methods to solve (11.14). 
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u(+1,t) = 0 (11.16) 
u(x,0) = £(x), au (x,0) = g(x) (11.17) 


is well posed. Legendre polynomial solution of (11.15-17) is 
semi-bounded and, hence stable (see Sec. 7). However, we have not 
yet been able to prove that Chebyshev solution of this problem is 
ever algebraically stable. The techniques of Sec. 8 prove that 


if the boundary conditions (11.16) are replaced by the characteristic 


conditions 


du(-l,t) , du(-1,t) = 9g, du(l,t) _ du(l,t) = 9, 
ot ot at at 

the scheme is algebraically stable. However, we have not yet 
been able to prove this result for (11.16). However, it is 
reassuring to note that we have solved the Cheybshev-spectral 
approximations to (11.15-17) and find no evidence of lack of 
convergence. Indeed, the Chebyshev methods work just as well as 
they do for (11.1-2). Thus, it is not the spectral methods that 
run into difficulty on higher-order equations, but just our 


methods of analysis. 
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12. Advective-Diffusion Equation 


In this section, we consider spectral methods for the 


advective-diffusion (‘linearized Burgers') equation 


2 


Qu(x,t) 4 yduliet) © 24 4 ¢€(x,¢) ta <1) (12.1) 
ot ox axe a = 

u(-1,t) = 0, u(l,t) = 0 (12.2) 

u(x,0) = g(x) (12.3) 


Eq. (12.1) is parabolic so boundary conditions should be applied 

at both x = -l and x=+#+l. When v is small, the boundary 
condition applied at x = +l (assuming U > 0) has an interesting 
effect on the stability of the spectral methods. 

To begin, we remark that the analyses of Sects. 7-8 can be 
extended to show that, as N+, N-term Legendre and Chebyshev 
approximations to (12.1-3) are stable and convergent. 

For example, Chebyshev-Galerkin approximation is stable 


because (12.1-2) and (7.3) imply that 


1 2 : 2 1 2 
a | - dx < | u dx - v | a ax 
he we) ees er (1=x2) 372 “1 (1-x2) 272 
10] 1-v/U uw? P 
rie U | . “Gr a. ore x 
a - /1-v/0 (1-x") 
2 1 2 
< = | ¥ dx (12.4) 
~1 A-x? 


so the approximation is semi-bounded. 


-235- 


a — 


However, for finite N, there may be difficulty integrating 
the resulting spectral equations. With Legendre polynomials, 


* 
Galerkin approximation Ly to (12.1-3) satisfies Ly +L. < 0 


N 
so there is no difficulty with time integrations (although the 
solution may not be accurate unless N_ is large enough. 
On the other hand, Chebyshev-spectral solution of (12.1-3) 
encounters the following curious behavior when v is small. 
If v/U is small and N_ is not too large, the Chebyshev~spectral 
approximations Ly to (12.1-3) have eigenvalues with positive 


real parts. In Table 12.1, we list values of N for various 


crit 


values of v/U; for N<N Ly for Chebyshev-tau 


eric ‘ 
approximation to (12.1-3) has eigenvalues with positive real 
parts. Since these eigenvalues may have moderately large real 
Parts [they can be as large as u2/2v by (12.4)], there may 
be rapid growth of errors and numerical solution of the 
Chebyshev~-spectral equations may appear unstable and divergent. 


For N>N there are no eigenvalues of Ly with positive 


erie. * 
real parts so the spectral equations are stable. 

The origin of this temporal instability is the outflow 
boundary layer at x = +l ; when U > 0, the solution to 
(12.1-3) develops a region of rapid change of width roughly 


v/U near x= +1 as t increases. Since roughly 3(u/v) 172 


Chebyshev polynomials are required to resolve a boundary layer 


of width v/U [see (3.50)], we expect that Noi, # 3 (v/v) 27? 
so UN 3/0 + a In fact, as shown in Table 12.1, the 


2 


ecie’ ¢4. [Since the Chebyshev norm 


1/4 


criterion is actually vN 


of exp(-Utd/3x) is roughly N (see Sec. 8), we expect 


that the proper scaling of N 


7/4 
' vn 2/5 70 $ 1.3. As shown in Table 12.1, this modified scaling is 


orae. «4 better represented as 


more nearly satisfied for the range of Vv considered. ] 
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TABLE 12.1 


ia ca 


| if - 
2 7/4 
v/U Norit Worit/? Worit/U | 
1.0 x 107? 15 2.25 1.14 
« | | 
| ' 2.5% 107° 35 3.06 1.26 
1.0 x 107? 61 3.72 1.33 
Q 
/ 6.0 x 1074 81 3.94 er 
| 
-4 
| 4.0 x 10 101 4.08 1.29 
a 2 eee ER HR eee RO ee en fe 
- Table 12.1 Critical values Norit of the number of Chebyshev 
polynomials necessary that the tau approximation to the operator 
-Udu/ox + va7u/ax with u(+l) = 0 have no eigenvalue with 
c 


positive real parts. Also listed are the inverse 'grid Reynolds 


7/4 /v 


/U and the parameter vNi i; 


2 
' 
number VWNorit 


If Chebyshev-spectral approximations to (12.1-3) are 


solved using fractional time-step methods, the temporal 


instability for N< N appears in a unique way. Define 


crit 
the operator Ay as an N-mode Chebyshev approximation to the 
operator -Udu/3x with the boundary condition u(-1) = 0 

and the operator By, as an N-mode Chebyshev approximation 

to the operator va7u/ax" with u(+l) = 0. Then the 
evolution operator of (12.1-2) is exp [ (Ay+B,) t] soa 


fractional step method involves the splitting 


duy/dt = 9 Uy/ ot + doy /dt where 
9 Uy/ot = AyYy ’ 9 2Uy/ dt = ByYy ‘ 


For any values of v and U> 0, the fractional step 

,U,/3t is algebraically stable since ||exp A,t|| = o(nt/4) 

(see Sec. 8), while the fractional step 9oUy/dt is stable 

since ||exp Byt || <2 (see Sec. 7). Nevertheless, 

l|exp [ (Ay+B,.) t] || can grow rapidly with t. The reason is 
N 

Jexp[(Ay+B,)t] |] < |lexp Aytl| llexp Bytll — - 


The Lie formula (5.8) does ensure that 


that Ay and B do not commute so it is not true that 


lJexp((Ay+By) t] || < — lJexp(Ayt/n) ||" llexp (Byt/n||", 


renee ct a ore ee oan mens reenact aeteemnccensmanceiee 4h 


However, as n+, 
llexp(A,t/n) || - 1 ~ CN7t/n 
with C> 0 (see Sec. 8) so 
llexp (Ayt/n) ||" ~ exp(CN7t) >> 21 (n+ ©) ; 


Therefore the Lie formula gives only the very weak upper bound 
2 
l|exp [ (A\+B,) t) || < exp(CN*t) . 


In summary, Chebyshev-spectral approximations to (12.1-3) 
give fractional step methods such that each fractional step 
is algebraically stable while the total step is unstable 


unless N>N 


crit 
If the boundary conditions (12.2) are replaced by 


“eRe “= @, au (41,t)) = y (12.4) 


when U > 0, the criterion for temporal stability is relaxed 


significantly. As shown in Table 12.2, the value of vn? /U 


crit 
is decreased to roughly 1.6. However, the growing modes that 


appear when N < N are much tamer than those that appear 


crit 
when the boundary condition u(+l,t) = 0 is applied, so 
accurate time integrations can still be obtained when 


vn?/u % 0.01 (see Haidvogel 1977). 
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TABLE 12.2 


Table 12.2. Critical values N of the number of 


erit 
Chebyshev polynomials necessary that the tau approxima- 


tion to the operator -Udu/dx + va7u/ax? with 
u(-1) = 0, du(+l)/ax = 0 and U>O have no 


eigenvalues with positive real parts. The parameter 


7/4 


VWNorit 


/U is also listed. 
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13. Models of Incompressible Fluid Dynamics 


The Stokes equations for low Reynolds number, two- 


dimensional incompressible flow are 


@ 
<+ 


Q' 
(+ 


Vev=0, 
where v is the velocity field, p is the pressure, and v is 
the kinematic viscosity. With the boundary conditions that 
v= 00n rigid stationary boundaries, the problem (13.1) is 
well posed for any v > 0. An equivalent formulation is given 


by the vorticity-streamfunction equations 


a vv7z ’ 


2 
Go = Nev» 


(13.2) 


obtained by taking the curl of the Stokes equations (13.1). Here 
y is the streamfunction defined by V = (-3/dy,3p/ax) and ¢ 
is the vorticity. 


A one-dimensional model of (13.2) is 


SE my 2 = 
= ae (<1 < x22, t 990), (13.3) 
2 
=. oS + “ 
. = <a el < eS Des (13.4) 


On stationary rigid walls, the boundary conditions for (13.3~-4) are 
w(x,t) = Vy (x,t) =z 0 (x = +1). (13.5) 


There is one subtlety in the application of spectral methods 
to (13.3-5) that does not appear directly when the primitive 
equations (13.1) are used. It is necessary to use some care 


to avoid unconditional numerical instability with the Chebyshev- 


tau method. 


The most obvious way to use the tau method to solve 


(13.3-5) is to substitute (13.4) into (13.3) and solve 


v t“L< =< 3, ¢ > 0) (13.6) 


xxt ~ YY¥xxxx 


by expanding y(x,t) in the Chebyshev series 
N 
v(x,t) = ) ai (t)T (x) . (13.7) 
n=0 


Denoting by a) the Chebyshev expansion coefficients of aty/axt 


(see A.20), the tau equations for (13.5-6) are 


da {?) 

n (4) 
“Te = a, (0 S n < N-4, <> 0), (13.8) 
N N 

Pipa, = EsGiy ns, = 8: (13.9) 
n=0 n=0 


Unfortunately, this method for solution of (13.3-5) is 
unconditionally unstable as N+. [In Table 13.1, we list 


the largest positive eigenvalue \ eae of (13.8-9); there is 


a solution of (13.8-9) that grows like a, (t) = a, (O)exp (A 
4 


max*) ° 


Since i grows like N° as N > ~, errors also grow rapidly 


max 
as N + ~ for fixed t. This method is unusable for time- 


dependent calculations. 


In Table 13.1, we also list the values of An for fn = 2,5) 
where the eigenvalues of (13.8-9) are ordered according to 
14! < |ro | < se. . The exact eigenvalues of (13.3-5) are 
found by seeking solutions of these equations of the form 
w(x,t) = w(x)exp(At), t(x,t) = t(x)exp(At). It may be easily 
verified that the exact eigenvalues of (13.3-5) are given by 

2 


A = <u" with wp = nt or uy any nonzero solution of the transcendental 


equation tan vy = up. The exact values of Ay and As are also listed 
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Table 13.1 


om 4s re 
c 10 -9.8696598 -189.63800 
15 -9.8696044 - 89.54550 
/ 20 -9.8696044 - 88.86244 
| -9.8696044 - 88.86244 
© | 30 -9.8696044 - 88.86244 
35 -9.8696044 - 88.86244 
40 -9.8696044 - 88.86244 
Exact -9.8696044 88.86244 


4,272. 
29,439. 
111,226. 
294,697. 
652,722. 
1,255,298. 
2,215,880. 


Table 13.1. Eigenvalues of the tau approximation (13.8-9) 
to (13.6-7). The N-4 eigenvalues are ordered so that 


last < [Agl & eee s [Ayag | -ALL the eigenvalues are real. 
The largest positive eigenvalue Min ae = i 


AE tee atl Mee 


in Table 13.1. Evidently, even though (13.8-9) is unstable as 
N > ~, it does a good job of reproducing the low-n modes; 
approximately yal Chebyshev polynomials are required to 
resolve the mode with eigenvalue ae Thus, this version 

of the tau method may be useful for eigenvalue calculations 
even though it is unconditionally unstable for the initial- 
value problem (13.3-5) (as evidenced by the spurious unstable 


modes with eigenvalues as large as dnax) * 


The tau method behaves similarly when applied to more 
difficult problems, like the Orr-Sommerfeld equation for 
linear stability analysis of incompressible plane-parallel 
shear flows. Low modes are given accurately by the analog 
of (13.8-9) (see Orszag 1971 ), but there appear spurious 
unstable modes with large growth rates. Similar spurious . 
unstable modes appear in finite~difference solution of the 


Orr-Sommerfeld equation (see Gary & Helgason 1970). 


There is a simple method to avoid the spurious unstable 
modes encountered by (13.8-9). The technique to be described 
below also eliminates the spurious unstable modes encountered 
in solution of the Orr-Sommerfeld equation. The idea is 
simply not to combine (13.3-4) into (13.6). Rather, we 


expand (x,t) as in 


N 
u(x,t) = qeghn (tn (13.10) 
and solve 
ab 
—p = vbl??) (oO <n <N-2), (13.11) 
= (2) < - 
bi, an (0 <n < N-2), (13.12) 


in addition to (13.9). Here we have dropped two equations 
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t) 


(! 


from the Chebyshev modal equations that result from (13.3-4). 

The logic of this modification of the tau method is as follows. 
Application of (13.8) for 0 < n < N-4 is equivalent to 
application of (13.12) for 0 < n < N together with (13.11) for 

0 <n < N-4,. On physical grounds, we may expect that this 
procedure will lead to instability because the boundary conditions 
y = 0 at x = +1 should be imposed on (13.4) not (13.3), while 


the boundary conditions y. = 0 at x = +1 should be imposed on 


x 
(13.3) only when v > 0. On the other hand, when the system 

is truncated as in (13.11-12), each of the dynamical equations 
can play their proper role in adjusting the boundary conditions: 
the boundary conditions yp = 0 are imposed on (13.12) while 


the boundary conditions vo * 0 are imposed on (13.11). 


We shall now prove that (13.11-12) is stable for the 
special case in which N is even with Qont+1 * bong = 9 for 
all n, t > 0. In this case, y(x,t) and ¢(x,t) are even functions 


of x. To begin , we observe that (13.11) is equivalent to 


2 t 
tev 23 + bay (x) (-L<x<l,t>0) , 
while (13.12) is equivalent to 
2 
cix,t) = 29 + bomix) . 
ax 
Therefore, 
he 34 
yCUF ov + @.1." . 
at axe ax we 


Since ) is an even function of x, it follows by integration with 


respect to x that 
2 3 


wot = 3 tat. 13.13) 
axot v ay™y (13. 
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Also, since y(x,t) is a polynomial of degree N that satisfies 
Vy (*1,t) = 0, integration by parts gives 


1 


i] - 1 - 
be, VogTy (1-x?) tax = - ee xy,-/ (1=x?) } By (1-x?) “Sax = 0 


since ¥. and x, / (1-*") are polynomials of degree N-2 so they must 
be orthogonal to Ty (x). Therefore, taking the Chebyshev inner 


product of (13.13) and v(x, t), we obtain 


1 1 
-). > 
x f , v2 (1-x”) Tax = 2v J Ya ocxx Max < 0, (13.14) 


where the last inequality is established using the inequality 


derived in Example 7.1(v): 
1 
fun. (1-x?)“#ax <0 
-. ** = 
if u(x) i: a polynomial of degree N satisfying u(+l) = 0. The 


energy bound (13.14) proves stability of the tau approximation 


(13.11-12). 


Finally, let us discuss methods for the solution of 
the primitive equations (13.1) using Chebyshev tau approximations. 
A one-dimensional model that embodies the essential features 
of (13.1) is obtained by solving (13.1) within the slab -1 <x <1, 
-© < y < ~, with an assumed solution of the form 


y, vix,tje-"’), p= p(x,t)er*Y 


v= (u(x,t)e 
for some real wavenumber k. Let the Chebyshev expansion 
coefficients of u(x,t), v(x,t), p(x,t) be denoted as u(t), 
Vv, ft) Py, (t) (0 < n < N), respectively. Then an unconditionally 
stable, implicit fractional step method for the solution of (13.1) 


with a forcing term (£ (x, t)et*Y g(x, t)et*¥) added to the right side is 
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- et 2 
U, = u(t) + Att-pov + Fi(t)) (0 < m <¢ N-2), (13.15) 
V, = V(t) + dt[-ikp, + gi (t)] (O0<n<N), (13.16) 
a, (1) + ikv, = 0 (0<n<N), (13.17) 
Pa, = 2 cum 
us = (-1)"u. = 0 iO <n<N), (13.18) 
n=0 ” n=0 ” walt at 
u (ttre) - vatul?) (teat) = 0, (0 <n < N-2), (13.19) 
v(ttat) - vatv(?) (teat) = 0, (0 <n < NH2), (13.20) 
8 n 4 n 
Y ( l)vu(t+at) = J (1) vi (t+at) = 0. (13.21) 
n=0 n=0 
Here we use the notation that, for example, a represents 


the Chebyshev coefficients of uy, (Xrt)- The scheme (13.15-21) 
is based on backwards Euler time differencing; it is straight- 
forward to generalize (13.15-21) to other more accurate time 
differencing methods. 

The fractional step (13.15~18) involves computation 
of the pressure field by imposition of the incompressibility 
condition (13.17). Only the boundary conditions u(+l1,t) = 0 
are applied because this part of the time step is effectively 
inviscid so only the normal flow can be specified at the 
boundary. Thus, we drop (13.15) for n = N-1,N in favor of 
the two boundary conditions (13.18). The fractional step 
(13.19-21) involves the viscous term in (13.1) so boundary 
conditions are applied on both the normal velocity component 
u and the tangential velocity component v. Accordingly, the 
tau method involves dropping (13.19-20) for n = N-1,N in favor 


cf these boundary conditions. 
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eta " a — 


The system (13.15-21) is solved as follows. 


Multiplying 
(13.15) by ik and subtracting the result from the Chebyshev 


x-derivative of (13.16) gives 


v,()) ~ aki, = vit) ce - iku(t) + at{g, () (t) - ike, (t)] 


Substituting Vn = ia from (13.17) gives 


2- 


= (2) 
un -k UL 


a ul?) (t)-k7u (t) + dt [-ikg ‘?) (t) - x7£ (t)] 


(0 <n <N-2). (13.22) 


Eq. (13.22) with the boundary conditions (13.18) is of the 
same form as (13.19-20) with boundary conditions (13.21). 
These equations are best solved by the algorithm discussed 
at the end of Sec. 10. 


The stability analysis of (13.15-21) is as follows. 


with fn =g.=0 for all n. Therefore, the solution of (13.22) 


The evolution of a perturbation is governed by (13.15-21) 
* 


is uy = u(t) for all n. Also, vn = v(t) for all n, Finally, 
the implicit scheme (13.19-21) is an unconditionally stable 
scheme for solution of the heat equation. This proves that 


(13.15=-21) is unconditionally stable. 


The methods discussed in this section extend to give 
stable methods for solution of the nonlinear Navier-Stokes equations. 
For example, if the forcing term (f,g) in (13.15-16) is chosen 
to be the nonlinear terms of the Navier-Stokes equations, our 


analysis shows that stability of (13.15-21) is determined by stability 


restrictions on the nonlinear terms alone. 


14. Miscellaneous Applications of Spectral Methods 


In this Section, we survey some special topics regarding 


spectral methods. Some of these topics are still under active 


investigation, so the results reported here are very incomplete. 


Complicated Geometries 

There are two ways that spectral methods can be used 
to solve problems in complicated geometries without introducing 
basis functions that are special to the geometry and, therefore, 
unwieldy and inefficient to use. The two methods are mapping 
and patching. 

Mapping involves transforming the complicated domain 
into a simpler one by means of a coordinate transformation. 
Spectral methods are then applied in the simple geometry 
using the techniques discussed in earlier sections. For 
example, if we wish to solve the heat equation 


« u(x,y,t) = V2u(x,y,t) (14.1) 


in the two-dimensional domain 
“l<x<1, -f£(x) sy < £(x) 
for some given function f(x) with the boundary conditions 
that u = 0 on the boundary of the domain, we would proceed 
as follows. First, we make the coordinate transformation 
2 = y/f (x) fw). $ 8S. 4) (14.2) 
and rewrite (14.1) as 


x u(x,z,t) = (se - e so) 7ulx,2,t) + f —z u(x,2,t) 


3 
fhe zed, 2 <8 Ss 4) (14.3) 


Then, we expand u(x,z,t) in a double Chebyshev series 
and integrate (14.3). For this purpose, a hybrid numerical 
scheme is suggested in which time differencing is stabilized 


by a semi-implicit method (see Sec. 10) in which a simple 


q diffusion operator is added and subtracted from (14.3). The 
simple diffusion operator is then evaluated using a tau 
j method (because the tau method is simplest when no complicated 


nonlinearities or nonconstant coefficient terms are involved); 


the remaining nonconstant coefficient term in (14.3) is then 
evaluated using fast Fourier transforms and the collocation 
method. The result is an efficient and accurate method for 
solution of (14.1). 


Techniques like those just described have been applied 


at a variety of problems with much success. If a convenient 
coordinate transformation is available, the mapping technique 
combined with appropriate spectral methods may be expected 
to be very useful. 

The idea of patching is that if the geometry is the 
union of several simpler geometries (like an L-shaped region) 
then spectral approximations can be formulated in each of 
the simpler domains and then patched across the boundaries by 
requiring that the solution (and an appropriate number of 
derivatives) be smooth. When this technique is applied 
together with the mapping technique discussed above, it is 
possible to devise spectral shock-fitting methods for the 
solution of compressible flow problems. These methods require 
much further investigation to judge their usefulness in practical 


problems. 
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Poisson's Equation in Two and Higher Dimensions 


The Chebyshev tau equations for Poisson's equation 


vu = f in the square -l<x<l, -l<y<l are 


u (270) 4 ylOr2) 2 £  (OsnsN-2,0<mmM-2), (14.4) 


while the Dirichlet boundary conditions u = 0 are 


N n 

25 (th) Wim = 9 (O<m<M) , (14.5) 
N 

wig th nm = 0 (O<n<N) . (14.6) 


Here we expand u(x,y) and f(x,y) in the double Chebyshev series 
N M u 
u (x,y) t nm 


and we denote the Chebyshev expansion coefficients of gP*4y /9xPayF 


by * alll The 2N+2M+4 boundary conditions are not all linearly 
independent; there exist four linear relations among them, namely 


N M 
n m - 
je Pie (+1) oy 0. (14.8) 


Thus, (14.4-6) gives (N+1) (M+l1) equations for the (N+1) (M+1) 


unknowns u_,, (0<n<N, O<m<M). 


Using (10.7) [or (A.20)], the system (14.4-6) can be reduced 
to a block tridiagonal matrix equation modified by extra full 
rows corresponding to the boundary conditions (14.5-6). These 
equations can be solved by standard block tridiagonal algorithms 
in order nom or order nM? operations. If Poisson's equation must 
be solved several times with the same values of N and M but different 


functions f(x,y), it is more efficient to apply alternative methods. 
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A method to solve (14.4-6) in order n2m operations (with 
a preprocessing stage that requires order n> operations) is 


as follows. First, we find the N-2 eigenvalues A and eigen- 


vectors ®np (p = 0,...,N-2) of the equations 
(zy sa 
enp = r5enp (0 <n < N-2) 
N 
n 
+1 = 0. 
a "enp 


The eigenvalues Ae are all negative as proved in Example 7.3(ii) . 


Then we form the (N+l)x(N+l) matrix E whose elements are 


aa < N- 
Enp enp (0 <n <N, 0 < p < N=2) 
a ee, iL 

= <n< 
En ,N Sh, eee SB) 


and compute the inverse matrix D = E °. Since the boundary 


conditions (14.5) are satisfied by Unm! it follows that 


> 
u = ev (14.9) 
nm p=0 np pm 
for suitable Vom for all n,m. Therefore, setting 
N 
i. * te (D) ontnm (OSPSN-2, O<MSM-2), (14.10) 
it follows that (14.4-6) become 
(0,2) _ <p <N- <m<M- 
AoYpm * Ypm’ = Spm (OSPSN-2, O<m<M-2) (14.11) 
M m. 
k (41)"v__ = 0 (0<p<N-2). (14.12) 
m-0 = “| ; 


Eqs. (14.11-12) may be solved efficiently (in order NM operations) 
for Yom using the algorithm discussed at the end of Sec. 10. 

Once Yom is found, Uam 
total operation count is order n2M [from the two matrix multiplies 


may be reconstructed from (]4.9). The 
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(14.9-10)]. 
The solution of Poisson's equation by the Chebyshev 
series method outlined above is very competitive with finite- 
difference solution using fast Poisson solvers. Zang & 
Haidvogel (1977) present a number of comparisons of the Chebyshev 


methods and fast Poisson solvers. 


There are two further complications that may arise in 
elliptic boundary-value problems. First, the elliptic equation 
may have nonconstant coefficients or may even be nonlinear. 

Here we recommend that spectral equations be solved using 

h relaxation methods of the kind advocated by Concus & Golub (1973), 

F in which the heart of the algorithm is the fast, efficient 
solution of Poisson-like equations. Second, the geometry 


may be more complicated than a box. In this case, we recommend 


Ww 


the implementation of capacitance matrix techniques (or 
equivalent Green's function techniques) in which the problem 


to be solved is imbedded in a simpler geometry, like a box 


« 
(see Buzbee et al 1971). Again, the heart of the algorithm 
is the fast solution of Poisson's equation using (14.9-12). 
C : ; Pact 
Coordinate Singularities 
When spectral methods are applied to problems in 
C cylindrical or spherical geometries, their formualtion may 
require special care at the coordinate singularities. These 
"pole problems' have been extensively investigated (Orszag 1974, 
C Tang 1977). As a simple example of these effects, let us 
consider the computation of the eigenvalues of Bessel's equation 
using the Chebyshev tau method (Metcalfe 1974). The problem is 
C -253- 


to find the eigenvalues and eigenfunctions y(x) of 


2 
« ere a) 
rg? a? Ay (14.13) 
subject to the conditions that y(1) = 0 and that y(x) be finite 
for 0 <x <1. The exact eigenvalues are related to the zeros of the 


a 
Bessel function Jy? » = JInp 


where Jy Gpp)=0- od a 
When n is even, the eigenfunctions of (14.13) are even 
functions of x; when n is odd, the eigenfunctions are odd. 
This fact suggests that we represent the solution to (14.13) 
in terms of series of even Chebyshev polynomials when n is even 


and odd polynomials when n is odd. Thus, for n odd we write 
M 


y(x) = ) VeTe th) (14.14) 
m=1 
In Table 14.1, we list numerical values for the smallest eigenvalue 
of (14.13) with n = 7 using the series (14.14), the boundary 
condition y(1) = 0, and the Chebyshev tau method. The convergence 
of this method, while very impressive as M increases, is slowed 
by the coordinate singularity of (14.13) at x = 0. In general, 


series of the form (14.14) behave like x as x +0. In this 


case the terms y'/x and y/x? are singular at x = 0. The true 


} 

} 

i 

| 

i 
eigenfunctions Jo (j,7*) behave like x! as x * 0, as may easily 
be shown using Frobenius’ method, so none of the terms of (14.13) 
are in fact singular for the exact eigenfunctions. 

it is possible to improve the convergence of (14.14) by 

imposing additional 'pole conditions', like y'(0) = 0. When 

| y'(0) = 0 in the series (14.14), the terms of (14.13) are 

| individually nonsingular. In Table 14.1, we also list numerical 

values of the smallest eigenvalue of (14.13) with n = 7 and 


| the two boundary conditions y(1) = 0, y'(0) = 0 applied. There 
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Table 14.1 


M Ay with y(1)=0 
ze, Pao 
14 169.111983340 
18 126.557832251 
22 122.991799598 
| 26 122.908250800 
Exact 122.907600204 


Table 14.1. Smallest eigenvalue of (14.13) with n = 7 
obtained using (14.14) and the Chebyshev tau method. 
M is the number of Chebyshev polynomials. 
boundary condition y'(0) = 0 is a pole constraint at 
the singular point x = 0 of (14.13). 


~124.001290649 


122.895944051 
122.907620295 
122.907600279 
122.907600204 
122.907600204 


Ay with y(1)=y' (0)=0 


The extra 


is clearly a dramatic improvement in the rate of convergence. 

It is also possible to make the problem less sensitive to 

pole properties near the origin by first multiplying (14.13) 

by x? to eliminate explicitly singular terms and then applying 
the tau method. The results of the latter trick are essentially 
the same as applying the pole condition y'(0) = 0 directly to 
(14.13). 

If pole conditions are not properly applied, it is possible 
to degrade significantly the accuracy of spectral computations. 
It is even possible to induce strong instabilities that are 
absent when proper pole conditions are applied. These matters 


are discussed in detail by Orszag (1974) and Tang (1977). 
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15. Survey of Spectral Methods and Applications 


In this Section, we give a brief survey of spectral 
' methods and some of their recent applications. There are 
five important features of spectral methods that should be 


i considered in their formulation and application. They are: 


i 
(i) Rate of convergence - If the solution to a problem 

is infinitely differentiable, then a properly designed 

: spectral method has the property that errors go to zero 

faster than any finite power of the number of retained modes. 

In contrast, finite-difference and finite-element methods 

yield finite-order rates of convergence. The important 
consequence is that spectral methods can achieve high accuracy 
with little more resolution than is required to achieve moderate 
accuracy. 

(ii) Efficiency - The development of fast transform 
methods permits spectral methods to be implemented with 
comparable efficiency to that of finite difference methods 
with the same number of independent degrees of freedom. 
However, since spectral methods typically require a factor 
of 2-5 fewer degrees of freedom in each space direction to 
achieve moderate accuracy (say, 5% error), the spectral 
computations can be considerably more effective. As the 
required accuracy increases, the attractiveness of spectral 
methods increases. 

(iii) Boundary conditions - As shown in earlier Sections 
of this monograph, the mathematical features of spectral 


methods follow very closely those of the partial differential 
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equation being solved. Thus, the boundary conditions imposed 
on spectral approximations are normally the same as those 
imposed on the differential equation. In contrast, finite- 
difference methods of higher order than the differential equation 
require additional ‘boundary conditions." Many of the 
complications of finite-order finite-diiference methods disappear 
with the infinite-order-accurate spectral methods. 

Another aspect of the treatment of boundary conditions 
by spectral tnethods is their high resolution of boundary 
layers. If the solution to a problem has a boundary layer 


of thickness ¢ , then only about 1/e% 


polynomials [see (3.50)] need 
be retained to achieve high accuracy. In contrast, finite- 
difference methods using equally spaced grid points would require 
about 1/e grid points to resolve such a boundary layer solution. 
Moreover, if a coordinate transformation is employed to improve 

the resolution of a boundary or internal layer of thickness ©€ , 

the errors of spectral methods are decreased faster than any 


finite power of ¢« as ©« >* 0. 


(iv) Discontinuities - Surprisingly, spectral methods 


do a better job of localizing errors than difference schemes 
and hence require considerably less local dissipation to smooth 
discontinuities. 

(v) Bootstrap estimation of accuracy - It is often 
possible to estimate the accuracy of spectral computations 
by examination of the shape of the spectrum. Thus, in computations 
of three-dimensional incompressible flows at high Reynolds numbers, 


the mean-square vorticity spectrum must not increase abruptly at 
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large wavenumbers (small scales); if the vorticity spectrum 
decreases smoothly to 0 as wavenumber increases, it is safe 
to infer that the calculation is accurate. On the other hand, 
similar criteria for finite-difference methods can be very 
misleading. 
Let us now survey some applications of spectral methods 
to incompressible fluid dynamics. We shall classify the method 


according to the boundary conditions and geometry. 


(i) Periodic boundary conditions in Cartesian coordinates - 


Here Fourier series are appropriate. Spectral methods have 
been regularly used in three dimensions with 32 x 32 x 32 
modes and in two dimensions with 128 x 128 modes to simulate 
homogeneous turbulence. Most operational codes now use 
pseudospectral (collocation) methods because aliasing errors 
are usually small. The key fast transform methods are described 
in detail by Orszag (1971c). 
More recently, more ambitious spectral codes have been 
developed. The KILOBOX code employs 1024 x 1024 Fourier modes 
in two dimensions while the CENTICUBE code uses up. to 
128 x 128 x 128 modes in three dimensions. These high resolution 
codes are now being used to study fundamental questions 
regarding high Reynolds number turbulence, including the structure 
of inertial ranges. 
7 (ii) Rigid boundary conditions in Cartesian coordinates - 
Here Chebyshev polynomials should be employed. Typical 
applications to date include numerical studies of turbulent 
shear flows and boundary layer transition. Pseudospectral 


methods are used, with Chebyshev polynomials particularly 
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convenient because fast Fourier transform methods can be 
applied. 

(iii) Rigid boundary conditions in cylindrical geometry - 
Here Chebyshev polynomials should be used in radius, Fourier 
series in angle, and either Fourier or Chebyshev series in 
the axial direction (depending on boundary conditions). Some 
technical aspects of the implementation of Chebyshev series 
in radius, including pole conditions, is discussed by Orszag 
(1974). Applications to date include studies of transition 
in circular Couette flow and pipe Poiseuille flow. In particular, 
it should be emphasized that Chebyshev polynomial expansions 
are much better suited for serious numerical work than the 
apparently more natural choice of Bessel function expansions 
in radius. There are two reasons: Chebyshev series converge 
faster to general functions regardless of their boundary 
conditions and Chebyshev-spectral methods can be implemented 
efficiently by fast transform methods. 

(iv) Problems in spherical geometry - Here surface 
harmonic expansions, generalized Fourier series, and ‘associated' 
Chebyshev expansions all have attractive features. A 
detailed discussion of these methods is outside the scope of 
this monograph, but roughly speaking generalized Fourier series 
permit the most efficient transform methods to be developed 
followed by associated Chebyshev expansions and then surface 
harmonic expansions but surface harmonic expansions are best 
with regard to the pole problem. A variety of applications 


of these methods to global atmospheric flows have been made. 


(v) Semi-infinite or infinite geometry - Here Chebyshev 
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expansions are best if the domain can be mapped or truncated 

to a finite domain without serious error. There are two cases 
here: additional boundary conditions may or may not be 

required at '‘infinity.' Here again the formulation of spectral 
methods thlicun closkiy the exact mathematics. If additional 
boundary conditions, like radiation or outflow boundary conditions, 
must be imposed on the truncated domain, then they should 

also be applied to the spectral method. On the other hand, 

if mapping without additional Weuneey conditions does not 


introduce a singularity in the exact equations, no boundary 


conditions at ‘infinity' are required in the spectral approximation. 
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Appendix. Properties of Chebyshev Polynomial Expansions 
The Chebyshev polynomial of degree n, Tf), is defined 


by 
T, (cosé) = cosnéa, (A.1) 


Thus, Tp (x) Se a Ty (x) = xX, To (x) = 2x*-1, T, (x) = 4x°~3x, 


Ty (x) = 8x4-8x+1, and so on. Some properties of Chebyshev 


polynomials are 
\ 
ia’ 
| T, (x) Lye | T! (x) |<n?, (A.2) 
Pp p-l 
<7, (42) = (#1)"*P “qn (n?-x?)/(2k+1), (A. 3) 
ax k=0 
aP 2p : 
lap T(x) |= O(n") (n=; p fixed), (A.4) 


n ee i 
Ty (0) = 0, TH (0) = (-1)"n. 
The following formulae relate the expansion coefficients 


an in the series 


f(x) = ano a, 7, (x) (|x| <1) 


to the expansion coefficients b, of 


Lf (x) = ce b, 1, (x) (|x| <2) 
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c for various linear operators L. We use the constants Cc 


| and d, defined by 


lo Co = 2, c= 9 (n<0), c=) (pd), 


d= 1 (m0), d= 0 (n<0). 


n 
Q 
Some formulae are: 
eo 
a) Lf = f£'(x): cb. = 2 pa (A. 6) 
ses ae P 
ptnodd 
Lf = £"(x): cb = J P(p°=n")a, (A.7) 
. p=n+2 
ptn even 
Li a +a,,) 
M Lf = xf(x): PY = 2'Sn-1°n-1 ‘ntl (A.8) 
2 oe {c a ~+(c.+c_,)a ta __.} 
LE = x f(x)t By |g fp 28n-2* nT on-1' Sn" Snt+2 (A.9) 
4 — "= 
Lf = x f(x)t bo = ae ee (A.10) 
iho ot eld et The he te,)a +(e, 17, ten 41 tens 2) neat ensa! 
I 
-n-1 
f(x)-f£(0) , es oo (A.11) 
oi. ammumine  Te  SM e oF 
p=ntl P 
ptn odd 
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S pon-2 
_ PeeetOeO Oe seb = 2 | kp-nd{-2) 2a (Aste 


al a p=n+2 P 
ptn even ; 
5 A.13) 
" (x) -£' (0 b = 4 } pa ( 
pe = Etre “nen p=nt+2 P ' 
p-n=2 mod 4 
t 
eo 
£' (x) -f£' (0) -£"(0)X , cb = 2{ ) (p-n+l)pa_ 
Lf = 2 m nn p=n+3 p 
- p-n=3 mod 4 
t 
‘ oo 
- } (p-n-1) pa, (A.14) 
p=n+5 
| P-n=l mod 4 
| 
= xf! : = es (A.15) 
L£ = xf'(x) : cb, = na, + 2 Bina p 
ptn even 
4 
| 
be = x2£* (x), b, = gl(n-Va,_yt (mtd) (14d, Cy) nen (A. 16) 
a2 
‘ 
+ 4 } pa } 
p=n+3 
ptn odd : 
; 
4 
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\ = « 2 2 
Lf = xf"(x): cb, = 2n(n+l)a,, + }  plp*-n -la, (4.17) 


p=n+3 
ptn odd 
re 5 2_,2 (x-T8) 
LE = x°£"(x): cnb, = n(n-lja, + y  plp*-n°-2)a (A.18) 
n Pp 
p=nt+2 
ptneven 
« £62) 
~ Lf ae 
-x 
ao 
ie with f£(#1)<0:c¢.b = -2 J] (p-nja (A.19) 
p=n+2 Pp 
ptneven 
Also, if we expand £'@) (x) as in 
aa . |, a) 
— f(x) = J a q @ (x) 
q , 
ax n=0 n n 
then 
(q) (q) 2 (q-1) .2 
Chn1 2 nel antl 2na, . (A. 20) 
ae as ee ca 
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